bianjieyuan.for

边界元简单模型计算 附有试验文档 地球物理边界元法数值模拟实验

      PROGRAM BEMIP
!	  PARAMETER (NS=25,NI=10,ND=35,M=45,NE=34,NS1=24)
	PARAMETER (NS=25,NI=12,ND=37,M=49,NE=36,NS1=24)
	  DOUBLE PRECISION FB(M,M),U0(M)
	  DIMENSION JN(2,NE),X(2,ND),W(3,4),SL(3,NE),ETS(NS1),
     *  FI(2),DI(2),F(ND,ND),D(ND,ND),U1(NS),U2(NS),US(NS),OM(ND),
     *  ROS(NS1),Q(2,4,NE),A(10),X0(NS1)
     	  DATA W/.930568,.069432,.173927,
     *         .669990,.330010,.326073,
     *         .330010,.669990,.326073,
     *         .069432,.930568,.173927/
       OPEN (4,FILE='BEM.TXT',STATUS='OLD')
	 READ(4,*) (A(I),I=1,NI)
	 READ(4,*) (X(1,I),I=1,ND)
	 READ(4,*) (X(2,I),I=1,ND)
	 READ(4,*) P1,P2,ET1,ET2
	 CLOSE(4)
	 I3=0
	 P11=P1
	 P22=P2
	 NQ=4
	 PI=3.1415926
	 CALL OMG(X,OM,NS,ND)
	 DO 600 I=NS+1,ND
	 I1=I-NS
600	 OM(I)=A(I1)
	 DO 20 L=1,NS-1
	 JN(1,L)=L
20	 JN(2,L)=L+1
	 DO 25 L=NS,NE-1
	 JN(1,L)=L+1
25	 JN(2,L)=L+2
	 JN(1,NE)=ND
	 JN(2,NE)=NS+1
	 CALL SLQ2(JN,X,W,SL,Q,ND,NE,NQ)
140	 DO 30 I=1,M
	 DO 30 J=1,M
30	 FB(I,J)=0
	 DO 35 I=1,ND
	 DO 35 J=1,ND
	 F(I,J)=0
35	 D(I,J)=0
	 DO 40 I=1,NS
	 DO 45 L=1,NS-1
	 CALL FIL(I,L,X,W,SL,Q,FI,ND,NE,NQ)
	 DO 45 N=1,2
45	 F(I,JN(N,L))=F(I,JN(N,L))+FI(N)
	 DO 50 L=NS,NE
	 CALL FIL(I,L,X,W,SL,Q,FI,ND,NE,NQ)
	 CALL DIL(I,L,X,W,SL,Q,DI,ND,NE,NQ)
	 DO 50 N=1,2
	 F(I,JN(N,L))=F(I,JN(N,L))+FI(N)
50	 D(I,JN(N,L))=D(I,JN(N,L))+DI(N) 
40	 CONTINUE
	 DO 55 I=1,NS
	 DO 60 J=1,NS
60	 FB(I,J)=-F(I,J)
	 FB(I,I)=FB(I,I)+OM(I)
	 DO 65 J=NS+1,NS+NI
65	 FB(I,J)=-F(I,J)
	 DO 70 J=NS+NI+1,NS+2*NI
	 J1=J-NI
70	 FB(I,J)=-D(I,J1)
55	 CONTINUE
	 DO 75 I=NS+1,NS+NI
	 DO 80 L=1,NS-1
	 CALL FIL(I,L,X,W,SL,Q,FI,ND,NE,NQ)
	 DO 80 N=1,2
80	 F(I,JN(N,L))=F(I,JN(N,L))+FI(N)
	 DO 85 L=NS,NE
	 CALL FIL(I,L,X,W,SL,Q,FI,ND,NE,NQ)
	 CALL DIL(I,L,X,W,SL,Q,DI,ND,NE,NQ)
	 DO 85 N=1,2
	 F(I,JN(N,L))=F(I,JN(N,L))+FI(N) 
85	 D(I,JN(N,L))=D(I,JN(N,L))+DI(N)
75	 CONTINUE
	 DO 90 I=NS+1,NS+NI
	 DO 95 J=1,NS
95	 FB(I,J)=-F(I,J)
	 DO 100 J=NS+1,NS+NI
100	 FB(I,J)=-F(I,J)
	 FB(I,I)=FB(I,I)+OM(I)
	 DO 105 J=NS+NI+1,NS+2*NI
	 J1=J-NI
105	 FB(I,J)=-D(I,J1)
90	 CONTINUE
	 DO 110 I=NS+NI+1,NS+2*NI
	 DO 115 J=NS+1,NS+NI
	 I1=I-NI
115	 FB(I,J)=F(I1,J)
	 J2=I-NI
	 FB(I,J2)=FB(I,J2)+2-OM(J2)
	 DO 120 J=NS+NI+1,NS+2*NI
	 J1=J-NI
	 I1=I-NI
120	 FB(I,J)=P22*D(I1,J1)/P11
110	 CONTINUE
	 DO 125 I=1,NS+NI
125	 U0(I)=-X(1,I)*P11
	 DO 130 I=NS+NI+1,NS+2*NI
130	 U0(I)=0
	 CALL SLSOE(FB,U0,M)
	WRITE(*,*)U0
	 IF(I3.EQ.1) GOTO 144
	 I3=1
	 DO 135 I=1,NS
135	 U1(I)=U0(I)
	 DO 136 I=1,NS-1
	 R=SQRT((X(1,NS-I+1)-X(1,NS-I))**2+(X(2,NS-I+1)-X(2,NS-I))**2)
136    ROS(I)=(U1(NS-I+1)-U1(NS-I))/R
	 P11=P1/(1.-ET1)
	 P22=P2/(1.-ET2)
	 GOTO 140
144    DO 145 I=1,NS
145    US(I)=U0(I)
       DO 150 I=1,NS
150    U2(I)=US(I)-U1(I)
	 DO 155 I=1,NS-1
	 X0(I)=(X(1,NS-I+1)+X(1,NS-I))/2.
155    ETS(I)=(U2(NS-I+1)-U2(NS-I))/(US(NS-I+1)-US(NS-I))
       OPEN(2,FILE='IP.TXT',STATUS='REPLACE')
	 WRITE(2,230)P1,P2,ET1,ET2
	 WRITE(2,260)X
	 WRITE(2,160)(U1(25-I+1),U2(25-I+1),US(25-I+1),I=1,25)
	 WRITE(2,165)(X0(I),ETS(I),ROS(I),I=1,24)
	 CLOSE(2)
230    FORMAT(5X,'P1=',F7.1,3X,'P2=',F7.1,3X,'ET1=',F7.3,3X,'ET2=',F7.3)
260    FORMAT(5X,'X=',F8.2,3X,'Y=',F8.2)
160    FORMAT(5X,'U1=',F12.5,2X,'U2=',F12.5,2X,'US=',F12.5)
165    FORMAT(5X,'X=',F7.2,3X,'ETS=',F12.5,3X,'ROS=',F12.5)
       STOP
	END
	
	
	SUBROUTINE SLQ2(JN,X,W,SL,Q,ND,NE,NQ)
	 DIMENSION JN(2,NE),X(2,ND),W(3,NQ),SL(3,NE),Q(2,NQ,NE)
	 DO 10 L=1,NE
	 J=JN(1,L)
	 K=JN(2,L)
	 SL(1,L)=X(2,K)-X(2,J)
	 SL(2,L)=-(X(1,K)-X(1,J))
	 SL(3,L)=SQRT(SL(1,L)**2+SL(2,L)**2)
	 DO 10 MQ=1,NQ
	 DO 10 N=1,2
10     Q(N,MQ,L)=X(N,J)*W(1,MQ)+X(N,K)*W(2,MQ)
       RETURN
	 END
	 
	 
	SUBROUTINE OMG(X,OM,NS,ND)
	 DIMENSION X(2,ND),OM(ND)
	 PI=3.1415926
	 DO 10 I=2,NS-1
	 XJ=X(1,I)-X(1,I-1)
	 YJ=X(2,I)-X(2,I-1)
	 XK=X(1,I+1)-X(1,I)
	 YK=X(2,I+1)-X(2,I)
10     OM(I)=1.-(ATAN(YK/XK)-ATAN(YJ/XJ))/PI
       OM(1)=1.
	 OM(NS)=1.
	 RETURN
	END 	 	    	 	 	 	 	     


	SUBROUTINE FIL(I,L,X,W,SL,Q,FI,ND,NE,NQ)
       DIMENSION X(2,ND),W(3,NQ),SL(3,NE),Q(2,NQ,NE),FI(2)
	 PI=3.1415926
	 FI(1)=0.
	 FI(2)=0.
	 XI=X(1,I)
	 YI=X(2,I)
	 XL=SL(1,L)
	 YL=SL(2,L)
	 DO 10 MQ=1,NQ
	   XQI=Q(1,MQ,L)-XI
	   YQI=Q(2,MQ,L)-YI
	   RQI=SQRT(XQI**2+YQI**2)
	   S=(XQI*XL+YQI*YL)*W(3,MQ)/(PI*RQI*RQI)
	   FI(1)=FI(1)+S*W(1,MQ)
10	   FI(2)=FI(2)+S*W(2,MQ)
	   RETURN
	 END 


	SUBROUTINE DIL(I,L,X,W,SL,Q,DI,ND,NE,NQ)
	 DIMENSION X(2,ND),W(3,NQ),Q(2,NQ,NE),DI(2),SL(3,NE)
	 PI=3.1415926
	 DI(1)=0.
	 DI(2)=0.
	 XI=X(1,I)
	 YI=X(2,I)
	 XYI=SL(3,L)
30	 DO 10 MQ=1,NQ
	 XQI=Q(1,MQ,L)-XI
	 YQI=Q(2,MQ,L)-YI
	 RQI=SQRT(XQI**2+YQI**2)
	 S=ALOG(1./RQI)*W(3,MQ)*XYI/PI
	 DI(1)=DI(1)+S*W(1,MQ)          
10     DI(2)=DI(2)+S*W(2,MQ)
40       RETURN
	 END 

	 SUBROUTINE SLSOE(A,B,N)
	 DOUBLE PRECISION A(N,N),B(N),C,D
	 N1=N-1
	 DO 100 K=1,N1
	 K1=K+1
	 C=A(K,K)
	 IF(DABS(C)-1.0D-15)1,1,3
1	 DO 7 J=K1,N
	 IF(DABS(A(J,K))-1.0D-15)7,7,5
5 	 DO 6 L=K,N
       C=A(K,L)
	 A(K,L)=A(J,L)
6      A(J,L)=C
       C=B(K)
	 B(K)=B(J)
	 B(J)=C
	 C=A(K,K) 
	 GOTO 3
7      CONTINUE
       D=0.
	 GOTO 300
3      C=A(K,K)
       DO 4 J=K1,N
4	 A(K,J)=A(K,J)/C
	 B(K)=B(K)/C
       DO 10 I=K1,N
	 C=A(I,K)
	 DO 9 J=K1,N
9      A(I,J)=A(I,J)-C*A(K,J)
10     B(I)=B(I)-C*B(K)
100    CONTINUE
       IF(DABS(A(N,N))-1.0D-15)11,11,101
11     WRITE(*,12)K
12     FORMAT('***SINGULARITY IN ROW',I5)
       D=0.
       GOTO 300
101    B(N)=B(N)/A(N,N)
	 DO 200 L=1,N1
	K=N-L
	K1=K+1
       DO 200 J=K1,N
200    B(K)=B(K)-A(K,J)*B(J)
       D=1.
	 DO 250 I=1,N
250    D=D*A(I,I)
300    RETURN
       END