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