PROGRAM LABPM1
C
C THIS PROGRAM SOLVES
2-D LAPLACE EQ. WITH SIMPLE BOUNDAY
C CONDITIONS USING
BOUNDARY POINT METHOD.
C THE FOLLOWING INPUT IS
REQUIRED:
C
C A0 = COEFFICIENT OF DIFFERENTIAL EQUATION.
C
C DD = DISTANCE BETWEEN THE CONTOUR OF SOURCES
AND THE
C BOUNDARY
OF THE BODY (ASSUMING HARMONIC CONTOUR)
C
C NB = NUMBER OF BOUNDARY POINTS (=NUMBER OF
SOURCE POINTS)
C
C XB,YB= X AND Y COORDINATES OF BOUNDARY POINTS. COORDINATES OF
C ONE
POINT PER LINE. N LINES.
C
C IBC= TYPE OF SPECIFIED BOUNDARY CONDITION AT
BOUNDARY POINTS
C (XB,YB).
IBC=1 FOR U AND IBC= 2 FOR Q.
C
C BC= VALUE OF BOUNDARY CONDITION SPECIFIED AT BOUNDARY
C POINTS(XB,YB)
ON THE BOUNDARY.
C
C NF= NUMBER OF POINTS IN THE DOMAIN
WHERE THE SOLUTION IS
C SOUGHT.
C
C XF,YF= X
AND Y COORDINATES OF DOMAIN POINTS. ENTER COORDINATES OF
C ONE
POINT PER LINE. NF LINES.
C
C NX,NY= X
AND Y COMPONENTS OF THE UNIT VECTOR NORMAL TO THE
C BOUNDARY
AT (XB,YB). NB LINES
C
PARAMETER(NB=50,NF=50)
DIMENSION
XB(NB),YB(NB),XS(NB),YS(NB),XF(NF),YF(NF)
DIMENSION
IBC(NB),BC(NB)
DIMENSION
A(NB,NB),B(NF,NB),C(NF,NB),D(NF,NB)
REAL NX(NB),NY(NB)
DIMENSION U(NF),Q1(NF),Q2(NF)
C
OPEN(1,FILE='LABPM1EX1.DAT')
OPEN(2,FILE='LABPM1EX1.OUT')
C
READ(1,*)A0
READ(1,*)DD
C
READ(1,*)(XB(I),YB(I),I=1,NB)
C
READ(1,*)(IBC(I),BC(I),I=1,NB)
C
READ(1,*)(XF(I),YF(I),I=1,NF)
C
READ(1,*)(NX(I),NY(I),I=1,NB)
C
XS=XB+DD*NX
YS=YB+DD*NY
C
WRITE(2,10)A0
WRITE(2,20)NB
DO I=1,NB
WRITE(2,30)XB(I),YB(I),IBC(I),BC(I)
END DO
C
CALL
MAT_A(A0,NB,XB,YB,NX,NY,XS,YS,IBC,A)
C
CALL
SOLVE(NB,A,BC)
C
CALL MAT_BCD(A0,NB,NF,XF,YF,XS,YS,B,C,D)
C
WRITE(2,40)NF
C
U=MATMUL(B,BC)
Q1=MATMUL(C,BC)
Q2=MATMUL(D,BC)
C
DO I=1,NF
WRITE(2,50)XF(I),YF(I),U(I),Q1(I),Q2(I)
END DO
C
10 FORMAT(///,'BOUNDARY POINT METHOD SOLUTION FOR TWO-DIMENSIONAL
+LAPLACE EQUATION',/,
'A0=',E12.4)
C
20 FORMAT(///,'BOUNDARY
CONDITIONS AT',I3,' BOUNDARY
POINTS',//,
+8X,'XB',12X,'YB',9X,'TYPE',5X,'VALUE')
C
30
FORMAT(2(2X,E12.5),6X,I1,2X,E12.4)
C
40 FORMAT(///,'SOLUTION AT ',
I2 ,' SELECTED POINTS',//,
+8X,'XF',12X,'YF',12X,'U',12X,'Q1',12X,'Q2')
C
50 FORMAT(5(2X,E12.4))
C
STOP
END
C
SUBROUTINE
MAT_A(A0,NB,XB,YB,NX,NY,XS,YS,IBC,A)
C
DIMENSION
XB(NB),YB(NB),XS(NB),YS(NB)
DIMENSION IBC(NB),A(NB,NB)
REAL NX(NB),NY(NB)
C
DO I=1,NB
DO J=1,NB
R2=(XB(I)-XS(J))*(XB(I)-XS(J))+(YB(I)-YS(J))*(YB(I)-YS(J))
IF(IBC(I).EQ.1)THEN
A(I,J)=0.5*LOG(R2)
C
ELSE
RX=XB(I)-XS(J)
RY=YB(I)-YS(J)
A(I,J)=A0*(RX*NX(I)+RY*NY(I))/R2
ENDIF
END DO
END
DO
C
RETURN
END
C
SUBROUTINE
MAT_BCD(A0,NB,NF,XF,YF,XS,YS,B,C,D)
C
DIMENSION
XF(NF),YF(NF),XS(NB),YS(NB),B(NF,NB),C(NF,NB),D(NF,NB)
C
DO I=1,NF
DO J=1,NB
R2=(XF(I)-XS(J))*(XF(I)-XS(J))+(YF(I)-YS(J))*(YF(I)-YS(J))
B(I,J)=0.5*LOG(R2)
C(I,J)=A0*(XF(I)-XS(J))/R2
D(I,J)=A0*(YF(I)-YS(J))/R2
END DO
END DO
RETURN
END
C
SUBROUTINE
SOLVE(N,A,RHS)
C
DIMENSION
A(N,N),RHS(N)
C
TOL=1.E-6
N1=N-1
DO 33 K=1,N1
K1=K+1
C=A(K,K)
IF(ABS(C)-TOL)23,23,25
23 DO 29 J=K1,N
IF(ABS((A(J,K)))-TOL)29,29,27
27 DO 28 L=K,N
C=A(K,L)
A(K,L)=A(J,L)
28 A(J,L)=C
C=RHS(K)
RHS(K)=RHS(J)
RHS(J)=C
C=A(K,K)
GO TO 25
29 CONTINUE
GO TO 30
25 C=A(K,K)
DO 26 J=K1,N
26 A(K,J)=A(K,J)/C
RHS(K)=RHS(K)/C
DO 32 I=K1,N
C=A(I,K)
DO 31 J=K1,N
31 A(I,J)=A(I,J)-C*A(K,J)
32 RHS(I)=RHS(I)-C*RHS(K)
33 CONTINUE
IF(ABS((A(N,N)))-TOL)30,30,34
34 RHS(N)=RHS(N)/A(N,N)
DO 35 L=1,N1
K=N-L
K1=K+1
DO 35 J=K1,N
35 RHS(K)=RHS(K)-A(K,J)*RHS(J)
GO TO 36
30 WRITE(2,24) K
24 FORMAT(//' ******
SINGULARITY IN ROW',I3,'******, BECAUSE OF:',/,
+'WRONG INPUT OF
COORD. & BOUNDARY COND. OR BAD CHOICE OF DD',//)
36 CONTINUE
C
RETURN
END
C