1000 DEFINT I-N:KEY OFF:OPTION BASE 1:GOSUB 1770
1010 PRINT"---  Bear & Verruijt - Groundwater Modeling"
1020 PRINT"---  Plane steady flow in triangular dam"
1030 PRINT"---  Finite Element Method":PRINT"---  Program 10.4"
1040 PRINT"---  Quadrangular elements, interactive input"
1050 PRINT"---  Graphical output added at Zurich, March 1987"
1060 PRINT"---  Non-homogeneity added at Zurich, March 1987":PRINT
1070 DIM X(240),Y(240),IP(240),F(240),U(240),V(240),W(240)
1080 DIM KP(240,10),P(240,10),NP(200,4),T(200):NZ=10
1090 DIM XJ(3),YJ(3),PJ(3,3),B(3),C(3),KS(4,3)
1100 DIM XB(21),XT(21),YT(21),CT(21),TR(20)
1110 INPUT"Width of dam at the bottom (m) ........ ";W
1120 INPUT"Inclination of slope left (degrees) ... ";AL
1130 INPUT"Inclination of slope right (degrees) .. ";AR
1140 INPUT"Water level left (m) .................. ";HL
1150 INPUT"Water level right (m) ................. ";HR
1160 INPUT"Number of columns of elements ......... ";NL
1170 IF NL>20 THEN NL=20:PRINT"   Maximum value : 20"
1180 INPUT"Number of rows of elements ............ ";NH
1190 IF NH>10 THEN NH=10:PRINT"   Maximum value : 10"
1200 GOSUB 1770:FOR I=1 TO NL
1210 PRINT"Permeability of column ";:PRINT USING "##";I;
1220 INPUT" ................ ";TR(I)
1230 NEXT I:PRINT:A=4*ATN(1)/180:AL=A*AL:AR=A*AR:E=HL/500
1240 N=(NH+1)*(NL+1):M=NH*NL
1250 CL=COS(AL)/SIN(AL):CR=COS(AR)/SIN(AR):IS=0
1260 FOR I=1 TO NL+1:A=(I-1)/NL:XB(I)=A*W:CT(I)=CL+A*(CR-CL)
1270 YT(I)=HL+.8*A*(HR-HL):XT(I)=XB(I)+CT(I)*YT(I):NEXT I
1280 FOR I=1 TO NL:K=(I-1)*(NH+1):L=(I-1)*NH:FOR J=1 TO NH:JA=L+J
1290 NP(JA,1)=K+J:NP(JA,2)=NP(JA,1)+NH+1:NP(JA,4)=NP(JA,1)+1
1300 NP(JA,3)=NP(JA,2)+1:T(JA)=TR(I):NEXT J,I
1310 PRINT"Calculation of pointer matrix":PRINT"   Element ..... ";
1320 FOR I=1 TO N:KP(I,1)=I:KP(I,NZ)=1:NEXT I:FOR J=1 TO M:PRINT J;
1330 FOR K=1 TO 4:KK=NP(J,K):FOR L=1 TO 4:LL=NP(J,L):IA=0
1340 FOR II=1 TO KP(KK,NZ):IF KP(KK,II)=LL THEN IA=1
1350 NEXT II:IF IA=0 THEN KB=KP(KK,NZ)+1:KP(KK,NZ)=KB:KP(KK,KB)=LL
1360 NEXT L,K,J:PRINT:LL=1
1370 FOR I=1 TO 4:FOR J=1 TO 3:K=I+J-1:IF K>4 THEN K=K-4
1380 KS(I,J)=K:NEXT J,I
1390 FOR J=1 TO NH+1:F(J)=HL:IP(J)=2:IP(N+1-J)=2:F(N+1-J)=HR:NEXT J
1400 IP(N)=0:FOR I=1 TO NL+1:K=(I-1)*NH+I
1410 FOR J=1 TO NH+1:L=K+J-1:Y(L)=(J-1)*YT(I)/NH
1420 IF IS=0 THEN F(L)=YT(I)
1430 X(L)=XB(I)+Y(L)*CT(I):NEXT J,I
1440 GOSUB 1790:LOCATE 1,1,0:IF LL<0 THEN 1630
1450 FOR I=N-NH TO N:F(I)=HR:IF Y(I)>HR THEN F(I)=Y(I)
1460 NEXT I:PRINT:PRINT"Generation of Matrix":PRINT"   Element ..... ";
1470 FOR I=1 TO N:FOR J=1 TO NZ:P(I,J)=0:NEXT J,I
1480 FOR J=1 TO M:PRINT J;:FOR KW=1 TO 4:FOR I=1 TO 3
1490 K=NP(J,KS(KW,I)):XJ(I)=X(K):YJ(I)=Y(K)
1500 NEXT I:B(1)=YJ(2)-YJ(3):B(2)=YJ(3)-YJ(1):B(3)=YJ(1)-YJ(2)
1510 C(1)=XJ(3)-XJ(2):C(2)=XJ(1)-XJ(3):C(3)=XJ(2)-XJ(1)
1520 D=ABS(XJ(1)*B(1)+XJ(2)*B(2)+XJ(3)*B(3)):IF D<EE THEN 1590
1530 D=T(J)/(4*D)
1540 FOR K=1 TO 3:KK=NP(J,KS(KW,K)):II=KP(KK,NZ):FOR LL=1 TO II:L=1
1550 KV=KS(KW,L):IF NP(J,KV)=KP(KK,LL) THEN 1570
1560 L=L+1:IF L<4 THEN 1550 ELSE 1580
1570 P(KK,LL)=P(KK,LL)+D*(B(K)*B(L)+C(K)*C(L))
1580 NEXT LL,K
1590 NEXT KW,J:PRINT:PRINT:GOSUB 1640
1600 LL=-1:FOR I=1 TO NL+1:K=I*(NH+1)
1610 A=F(K)-Y(K):YT(I)=YT(I)+A:IF ABS(A)>E THEN LL=1
1620 NEXT I:IS=IS+1:PRINT:GOTO 1400
1630 GOSUB 1880:SCREEN 0:PRINT:END
1640 PRINT"Solution of equations":PRINT"   Iteration ... ";:IT=1
1650 EE=.00001:EE=EE*EE:FOR I=1 TO N:U(I)=0:IF IP(I)>0 THEN 1670
1660 FOR J=1 TO KP(I,NZ):U(I)=U(I)-P(I,J)*F(KP(I,J)):NEXT J
1670 V(I)=U(I):NEXT I:UU=0:FOR I=1 TO N:UU=UU+U(I)*U(I):NEXT I
1680 PRINT IT;:FOR I=1 TO N:W(I)=0
1690 FOR J=1 TO KP(I,NZ):W(I)=W(I)+P(I,J)*V(KP(I,J)):NEXT J,I
1700 VW=0:FOR I=1 TO N:VW=VW+V(I)*W(I):NEXT I
1710 AA=UU/VW:FOR I=1 TO N:IF IP(I)>0 THEN 1730
1720 F(I)=F(I)+AA*V(I):U(I)=U(I)-AA*W(I)
1730 NEXT I:WW=0:FOR I=1 TO N:WW=WW+U(I)*U(I):NEXT I
1740 BB=WW/UU:FOR I=1 TO N:V(I)=U(I)+BB*V(I):NEXT I:UU=WW
1750 IT=IT+1:IF IT<=N AND UU>EE THEN 1680
1760 RETURN
1770 SCREEN 0:CLS:LOCATE 1,30,1:COLOR 0,7:PRINT" Finite Elements - 4 ";
1780 COLOR 7,0:PRINT:PRINT:RETURN
1790 SCREEN 2:CLS:SX=640:SY=190:SXY=2:XA=X(1):YA=Y(1):XB=XA:YB=YA
1800 FOR I=1 TO N:IF X(I)<XA THEN XA=X(I)
1810 IF X(I)>XB THEN XB=X(I)
1820 IF Y(I)<YA THEN YA=Y(I)
1830 IF Y(I)>YB THEN YB=Y(I)
1840 NEXT I:DX=XB-XA:DY=YB-YA:SX=SX/DX:SY=SXY*SY/DY:IF SY<SX THEN SX=SY
1850 SY=SX/SXY:FOR I=1 TO M:FOR J=1 TO 4:JJ=NP(I,J):K=J+1:IF K>4 THEN K=1
1860 KK=NP(I,K):IA=SX*X(JJ):IB=SX*X(KK):JA=190-SY*Y(JJ):JB=190-SY*Y(KK)
1870 LINE (IA,JA)-(IB,JB):NEXT J,I:RETURN
1880 A$=INPUT$(1):RETURN