PROGRAM INTEFG
C======================================================================
C                       BOUNDARY ELEMENT METHOD
C                              APPLIED TO
C              LAPLACE EQUATION (STEADY HEAT EQUATION )
C    THIS PROGRAM COMPUTES G AND F FUNCTIONS FOR A GIVEN SET OF DATA
C    PROGRAMMED BY EIJI FUKUMORI // SUNY AT BUFFALO  // 1984 SPRING
C======================================================================
      IMPLICIT REAL*8 ( A-H , O-Z )
      PARAMETER (MXE=60,INTEPT=4, ND=2)
      DIMENSION XE(ND),YE(ND),SAI(INTEPT),W(INTEPT),X(MXE),Y(MXE),
     * NODEX(MXE,ND)
C========================= INITIAL SET-UPS ============================
      C1 = - 1./ ( 8.* DATAN( 1.D0) )
      CALL GRULE ( INTEPT, SAI, W )
C========================= READING IN DATA ============================
      CALL INPUT (ND, MXE, MXI, NE, NODEX, X, Y )
C=============== COMPUTATION OF F AND G FUNCTIONS  ==============
C
      OPEN ( 1, FILE='FUNCTION.FG', STATUS='UNKNOWN' )
      WRITE (1,*) '  SOURCE-ELEMENT CURRENT-ELEMENT  F   G'
      DO ISOURCE = 1 , NE
      DO ICURREN = 1 , NE
      CALL MATRIX ( ND,MXE,C1,INTEPT,SAI,W,NE,NODEX,X,Y,XE,YE, 
     * ISOURCE, ICURREN )
      END DO
      END DO
      CLOSE (1)
      STOP 'NORMAL TERMINATION'
      END
C
C
      SUBROUTINE MATRIX (ND,MXE,C1,INTEPT,SAI,W,NE,NODEX, X,Y,XE,YE,
     *  ISOURCE, ICURREN  )
      IMPLICIT REAL*8 ( A-H , O-Z )
      DIMENSION X(MXE),Y(MXE), NODEX(MXE,ND),XE(ND),YE(ND),
     *  SAI(INTEPT),W(INTEPT) 
C .......     (XP,YP) = COORDINATES OF SOURCE POINT.
      XP = ( X(NODEX(ISOURCE,1)) + X(NODEX(ISOURCE,2)) ) / 2.
      YP = ( Y(NODEX(ISOURCE,1)) + Y(NODEX(ISOURCE,2)) ) / 2.
C ------ CURRENT ELEMENT
      DO  I = 1 , ND
      XE(I) = X(NODEX(ICURREN,I))
      YE(I) = Y(NODEX(ICURREN,I))
      END DO
C-------- INTEGRATION OF f AND G ON A GIVEN ELEMENT
      IF ( ISOURCE .EQ. ICURREN ) THEN
      CALL FINE ( ND, C1, XE, YE, GE, FE )
      ELSE
      CALL INTE ( INTEPT,ND, XP, YP, C1, XE, YE, SAI, W, GE, FE )
      END IF
      WRITE (1,*) ISOURCE, ICURREN , FE, GE
      RETURN
      END
C
C
      SUBROUTINE INPUT (ND,MXE,MXI,NE,NODEX, X,Y )
      IMPLICIT REAL*8 ( A-H , O-Z )
      DIMENSION X(MXE),Y(MXE),NODEX(MXE,ND)
      OPEN ( 1, FILE='BEMFG.DAT', STATUS='OLD' )
      READ (1,*) NE
      DO IEL = 1 , NE
      READ (1,*) I,(NODEX(I,J),J=1,ND)
      END DO
      DO I = 1 , NE
      READ (1,*) NODE, X(NODE), Y(NODE)
      END DO
      CLOSE (1)
      RETURN
      END
C
C
      SUBROUTINE FINE ( ND, C1, XE, YE, GE, FE )
      IMPLICIT REAL*8 ( A-H , O-Z )
      DIMENSION XE(ND), YE(ND)
      FE = 0.D0
      DX = XE(2) - XE(1)
      DY = YE(2) - YE(1)
      DS = DSQRT ( DX*DX + DY*DY )
      GE = DS*C1 * ( DLOG(DS/2.D0) - 1.D0 )
      RETURN
      END
C
C
      SUBROUTINE INTE (INTEPT, ND, XP, YP, C1, XE, YE, SAI,W, GE,FE )
      IMPLICIT REAL*8 ( A-H , O-Z )
      DIMENSION XE(ND), YE(ND), SAI(INTEPT), W(INTEPT)
      GE = 0.D0
      FE = 0.D0
      DX = XE(2) - XE(1)
      DY = YE(2) - YE(1)
      DS = DSQRT ( DX*DX + DY*DY )
C-------(RNX, RNY) REPRESENTS NORMAL VECTOR
      RNX =  DY/DS
      RNY = -DX/DS
C------- DETJ STANDS FOR DETERMINANT OF JACOBIAN MATORIX
      DETJ = DS/2.D0
      XM = ( XE(2) + XE(1) ) /2.D0
      YM = ( YE(2) + YE(1) ) /2.D0
C
      DO IGAUSS = 1 , INTEPT
      XGAUSS = DX/2.D0*SAI(IGAUSS) + XM
      YGAUSS = DY/2.D0*SAI(IGAUSS) + YM
      RX = XGAUSS - XP
      RY = YGAUSS - YP
      R = DSQRT ( RX*RX + RY*RY )
C-------- INTEGRATION OF G(R)
      GE = GE + DLOG(R) * W(IGAUSS)
C-------- INTEGRATION OF F(R)
      FE = FE + ( RX*RNX + RY*RNY ) / (R*R) * W(IGAUSS)
      END DO
      GE =  C1 * DETJ * GE
      FE = -C1 * DETJ * FE
      RETURN
      END
C
C
      SUBROUTINE GRULE ( N , SAI , W )
      IMPLICIT REAL*8 ( A-H , O-Z )
      DIMENSION SAI(N) , W(N)
      IF ( N .LT. 2 ) STOP'N<2'
      IF ( N .GT. 4 ) STOP'N>4'
      IF ( N .EQ. 2 ) THEN
      SAI(1) = DSQRT(3.D0)/3.D0
      W(1) = 1.D0
      SAI(2) = - SAI(1)
      W(2) = W(1)
      RETURN
      END IF
      IF ( N .EQ. 3 ) THEN
      SAI(1) = DSQRT(15.D0)/5.D0
      SAI(2) = 0.D0
      W(1) = 5.D0/ 9.D0
      W(2) = 8.D0/ 9.D0
      SAI(3) = - SAI(1)
      W(3) = W(1)
      RETURN
      END IF
      IF ( N .EQ. 4 ) THEN
      SAI(1) = 0.33998104358485D0
      SAI(2) = 0.86113631159405D0
        W(1) = 0.65214515486254D0
        W(2) = 0.34785484513745D0
      SAI(3) = -SAI(2)
      SAI(4) = -SAI(1)
      W(3) = W(2)
      W(4) = W(1)
      RETURN
      END IF
      RETURN
      END