1000 DEFINT I-N:KEY OFF:OPTION BASE 1:GOSUB 1760
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":PRINT
1050 DIM X(240),Y(240),IP(240),F(240),U(240),V(240),W(240)
1060 DIM KP(240,10),P(240,10),NP(200,4),T(200):NZ=10
1070 DIM XJ(3),YJ(3),PJ(3,3),B(3),C(3),KS(4,3)
1080 DIM XB(21),XT(21),YT(21),CT(21)
1090 INPUT"Width of dam at the bottom (m) ........ ";W
1100 INPUT"Inclination of slope left (degrees) ... ";AL
1110 INPUT"Inclination of slope right (degrees) .. ";AR
1120 INPUT"Water level left (m) .................. ";HL
1130 INPUT"Water level right (m) ................. ";HR
1140 INPUT"Number of columns of elements ......... ";NL
1150 IF NL>20 THEN NL=20:PRINT"   Maximum value : 20"
1160 INPUT"Number of rows of elements ............ ";NH
1170 IF NH>10 THEN NH=10:PRINT"   Maximum value : 10"
1180 PRINT:A=4*ATN(1)/180:AL=A*AL:AR=A*AR:E=HL/500
1190 N=(NH+1)*(NL+1):M=NH*NL
1200 CL=COS(AL)/SIN(AL):CR=COS(AR)/SIN(AR):IS=0
1210 FOR I=1 TO NL+1:A=(I-1)/NL:XB(I)=A*W:CT(I)=CL+A*(CR-CL)
1220 YT(I)=HL+.8*A*(HR-HL):XT(I)=XB(I)+CT(I)*YT(I):NEXT I
1230 FOR I=1 TO NL:K=(I-1)*(NH+1):L=(I-1)*NH:FOR J=1 TO NH:JA=L+J
1240 NP(JA,1)=K+J:NP(JA,2)=NP(JA,1)+NH+1:NP(JA,3)=NP(JA,1)+1
1250 NP(JA,4)=NP(JA,2)+1:T(JA)=1:NEXT J,I
1260 PRINT"Calculation of pointer matrix":PRINT"   Element ..... ";
1270 FOR I=1 TO N:KP(I,1)=I:KP(I,NZ)=1:NEXT I:FOR J=1 TO M:PRINT J;
1280 FOR K=1 TO 4:KK=NP(J,K):FOR L=1 TO 4:LL=NP(J,L):IA=0
1290 FOR II=1 TO KP(KK,NZ):IF KP(KK,II)=LL THEN IA=1
1300 NEXT II:IF IA=0 THEN KB=KP(KK,NZ)+1:KP(KK,NZ)=KB:KP(KK,KB)=LL
1310 NEXT L,K,J:PRINT:LL=1
1320 FOR I=1 TO 4:FOR J=1 TO 3:K=I+J-1:IF K>4 THEN K=K-4
1330 KS(I,J)=K:NEXT J,I
1340 FOR J=1 TO NH+1:F(J)=HL:IP(J)=2:IP(N+1-J)=2:F(N+1-J)=HR:NEXT J
1350 IP(N)=0:FOR I=1 TO NL+1:K=(I-1)*NH+I
1360 FOR J=1 TO NH+1:L=K+J-1:Y(L)=(J-1)*YT(I)/NH
1370 IF IS=0 THEN F(L)=YT(I)
1380 X(L)=XB(I)+Y(L)*CT(I):NEXT J,I
1390 GOSUB 1760:PRINT"Free surface after";IS;"cycle";
1400 IF IS<>1 THEN PRINT"s";
1410 PRINT:PRINT:FOR I=1 TO NL+1:K=I*(NH+1):A$="####.###"
1420 PRINT"  x = ";:PRINT USING A$;X(K);
1430 PRINT"  y = ";:PRINT USING A$;Y(K):NEXT I:IF LL<0 THEN 1620
1440 FOR I=N-NH TO N:F(I)=HR:IF Y(I)>HR THEN F(I)=Y(I)
1450 NEXT I:PRINT:PRINT"Generation of Matrix":PRINT"   Element ..... ";
1460 FOR I=1 TO N:FOR J=1 TO NZ:P(I,J)=0:NEXT J,I
1470 FOR J=1 TO M:PRINT J;:FOR KW=1 TO 4:FOR I=1 TO 3
1480 K=NP(J,KS(KW,I)):XJ(I)=X(K):YJ(I)=Y(K)
1490 NEXT I:B(1)=YJ(2)-YJ(3):B(2)=YJ(3)-YJ(1):B(3)=YJ(1)-YJ(2)
1500 C(1)=XJ(3)-XJ(2):C(2)=XJ(1)-XJ(3):C(3)=XJ(2)-XJ(1)
1510 D=ABS(XJ(1)*B(1)+XJ(2)*B(2)+XJ(3)*B(3)):IF D<EE THEN 1580
1520 D=T(J)/(4*D)
1530 FOR K=1 TO 3:KK=NP(J,KS(KW,K)):II=KP(KK,NZ):FOR LL=1 TO II:L=1
1540 KV=KS(KW,L):IF NP(J,KV)=KP(KK,LL) THEN 1560
1550 L=L+1:IF L<4 THEN 1540 ELSE 1570
1560 P(KK,LL)=P(KK,LL)+D*(B(K)*B(L)+C(K)*C(L))
1570 NEXT LL,K
1580 NEXT KW,J:PRINT:PRINT:GOSUB 1630
1590 LL=-1:FOR I=1 TO NL+1:K=I*(NH+1)
1600 A=F(K)-Y(K):YT(I)=YT(I)+A:IF ABS(A)>E THEN LL=1
1610 NEXT I:IS=IS+1:PRINT:GOTO 1350
1620 PRINT:END
1630 PRINT"Solution of equations":PRINT"   Iteration ... ";:IT=1
1640 EE=.00001:EE=EE*EE:FOR I=1 TO N:U(I)=0:IF IP(I)>0 THEN 1660
1650 FOR J=1 TO KP(I,NZ):U(I)=U(I)-P(I,J)*F(KP(I,J)):NEXT J
1660 V(I)=U(I):NEXT I:UU=0:FOR I=1 TO N:UU=UU+U(I)*U(I):NEXT I
1670 PRINT IT;:FOR I=1 TO N:W(I)=0
1680 FOR J=1 TO KP(I,NZ):W(I)=W(I)+P(I,J)*V(KP(I,J)):NEXT J,I
1690 VW=0:FOR I=1 TO N:VW=VW+V(I)*W(I):NEXT I
1700 AA=UU/VW:FOR I=1 TO N:IF IP(I)>0 THEN 1720
1710 F(I)=F(I)+AA*V(I):U(I)=U(I)-AA*W(I)
1720 NEXT I:WW=0:FOR I=1 TO N:WW=WW+U(I)*U(I):NEXT I
1730 BB=WW/UU:FOR I=1 TO N:V(I)=U(I)+BB*V(I):NEXT I:UU=WW
1740 IT=IT+1:IF IT<=N AND UU>EE THEN 1670
1750 RETURN
1760 CLS:LOCATE 1,30,1:COLOR 0,7:PRINT" Finite Elements - 4 ";
1770 COLOR 7,0:PRINT:PRINT:RETURN