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

C

C

C

C

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