PROGRAM BEM8CONQ
C======================================================================
CC                       BOUNDARY ELEMENT METHOD
CC                              APPLIED TO
CC              LAPLACE EQUATION (STEADY HEAT EQUATION )
CC                 ELEMENT TYPE: CONSTANT ELEMENT
CC    PROGRAMMED BY EIJI FUKUMORI // NAGOYA // MARCH 21, 2022
C======================================================================
      IMPLICIT REAL*8 ( A-H , O-Z )
      PARAMETER (MXE=60,MXI=20,INTEPT=4, ND=2)
      DIMENSION XE(ND),YE(ND),
     * SAI(INTEPT),W(INTEPT),
     * A(MXE,MXE),
     * NODEX(MXE,ND),
     * C(MXE),
     * X(MXE),Y(MXE),IELTYPE(MXE),BV(MXE),
     * QN(MXE),H(MXE),RHS(MXE),
     * XI(MXI),YI(MXI), HI(MXI), CI(MXI), QX(MXI), QY(MXI)
C========================= INITIAL SET-UPS ============================
C
      C1 = - 1./ ( 8.* DATAN( 1.D0) )
      CALL GRULE ( INTEPT, SAI, W )
C========================= READING IN DATA ============================
C
      CALL INPUT (ND,MXE,MXI,NE,NODEX,IELTYPE,BV,X,Y,NIPT,XI,YI)
CC      ...TYPE(I)=1 ---> TEMPERATURE PRESCRIBED
CC      ...TYPE(I)=2 ---> HEAT FLUX PRESCRIBED
C============ BOUNDARY VALUE ARRANGEMENT ==============================
C
      CALL BVARRANG ( MXE, NE, IELTYPE, BV, H, QN )
C=============== FORMATION OF MATRIX G,F  AND  VECTOR C  ==============
C=============== FORMATION OF MATRIX A  AND  VECTOR RHS  ==============
CC
      CALL MTXFORM (ND,MXE,C1,INTEPT,SAI,W,NE,NODEX,X,Y,XE,YE,
     *                   C ,IELTYPE,RHS,H,QN, A )
C===================== READY TO SOLVE  A . X  =  C  ===================
C
      CALL SYSTEM ( MXE , NE, A , RHS )
C====================== SORTING SOLUTION ==============================
C
      CALL SORTSOLN ( MXE, NE, IELTYPE, RHS, H, QN )
C================== INTERNAL POINTS   =================================
C
      CALL DOMAIN ( INTEPT,ND,MXE,MXI,C1,NIPT,NE,SAI,W,
     *            XI,YI, NODEX, X,Y, H,QN, HI, CI, XE, YE)
      CALL DOMAIN2 ( INTEPT,ND,MXE,MXI,C1,NIPT,NE,SAI,W,
     *            XI,YI, NODEX, X,Y, H,QN, QX, QY, XE, YE)

C================= PRINTING RESULTS   =================================
      CALL ECHOSOLN( MXE,MXI,NE,ND,NIPT,NODEX,X,Y,IELTYPE,H,QN,
     *               BV,XI,YI,HI,CI, C, QX, QY)
      STOP 'NORMAL TERMINATION'
      END
C
C
      SUBROUTINE SORTSOLN ( MXE, NE, IELTYPE, RHS, H, QN )
      IMPLICIT REAL*8 ( A-H , O-Z )
      DIMENSION IELTYPE(MXE),RHS(MXE),QN(MXE),H(MXE)
      DO  I = 1 , NE
      IF ( IELTYPE(I) .EQ. 1 ) QN(I) = RHS(I)
      IF ( IELTYPE(I) .EQ. 2 ) H(I) = RHS(I)
      END DO
      RETURN
      END
C
C
      SUBROUTINE BVARRANG ( MXE, NE, IELTYPE, BV, H, QN )
      IMPLICIT REAL*8 ( A-H , O-Z )
      DIMENSION IELTYPE(MXE),BV(MXE),QN(MXE),H(MXE)
      DO I = 1 , NE
      IF ( IELTYPE(I) .EQ. 1 ) H(I) = BV(I)
      IF ( IELTYPE(I) .EQ. 2 ) QN(I) = BV(I)
      END DO
      RETURN
      END
C
C
      SUBROUTINE MTXFORM (ND,MXE,C1,INTEPT,SAI,W,NE,NODEX,X,Y,XE,YE,
     *                   C ,IELTYPE,RHS,H,QN, A )
      IMPLICIT REAL*8 ( A-H , O-Z )
      DIMENSION C(MXE),X(MXE),Y(MXE),
     * XE(ND),YE(ND),SAI(INTEPT),W(INTEPT) , NODEX(MXE,ND),
     * IELTYPE(MXE),QN(MXE),H(MXE),RHS(MXE),A(MXE,MXE)
C=============== FORMATION OF MATRIX G,F  AND  VECTOR C  
C----------  CLEAR MATRIX G(I,J) AND F(I,J) AND C(I)
      DO  I = 1 , NE
      C(I) = 0.D0
      RHS(I) = 0.D0
      DO  J = 1 , NE
      A(I,J) = 0.D0
      END DO
      END DO
C .......     (XP,YP) = COORDINATES OF OBSERVATION POINT.
      DO IOBSERV = 1 , NE
      XP = ( X(NODEX(IOBSERV,1)) + X(NODEX(IOBSERV,2)) ) / 2.
      YP = ( Y(NODEX(IOBSERV,1)) + Y(NODEX(IOBSERV,2)) ) / 2.
      DO ICURREN = 1 , NE
      DO  I = 1 , ND
      XE(I) = X(NODEX(ICURREN,I))
      YE(I) = Y(NODEX(ICURREN,I))
      END DO
C------------------ INTEGRAL ON AN ELEMENT
      IF ( IOBSERV .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
C------------------ MATRIX FORMATION
C.............NON-FREE TERMS
      IF (IELTYPE(ICURREN) .EQ. 1 ) THEN
      A(IOBSERV,ICURREN) = GE
      RHS(IOBSERV) = RHS(IOBSERV) + FE * H(ICURREN)
      END IF
      IF (IELTYPE(ICURREN) .EQ. 2 ) THEN
      A(IOBSERV,ICURREN) = - FE
      RHS(IOBSERV) = RHS(IOBSERV) - GE * QN(ICURREN)
      END IF
      C(IOBSERV) = C(IOBSERV) + FE
      END DO
C.............FREE TERM
      IF (IELTYPE(IOBSERV) .EQ. 1 ) THEN
      RHS(IOBSERV) = RHS(IOBSERV) - C(IOBSERV) * H(IOBSERV)
      END IF
      IF (IELTYPE(IOBSERV) .EQ. 2 ) THEN
      A(IOBSERV,IOBSERV) = A(IOBSERV,IOBSERV) + C(IOBSERV)
      END IF
      END DO
      RETURN
      END
C
C
      SUBROUTINE INPUT (ND,MXE,MXI,NE,NODEX,IELTYPE,BV,X,Y,
     * NIPT,XI,YI)
      IMPLICIT REAL*8 ( A-H , O-Z )
      DIMENSION X(MXE),Y(MXE),IELTYPE(MXE),BV(MXE),NODEX(MXE,ND),
     * XI(MXI),YI(MXI)
      OPEN ( 1, FILE='BEM0.DAT', STATUS='OLD' )
      READ (1,*) NE
      DO IEL = 1 , NE
      READ (1,*) I,(NODEX(I,J),J=1,ND),IELTYPE(I), BV(I)
      END DO
      DO I = 1 , NE
      READ (1,*) NODE, X(NODE), Y(NODE)
      END DO
      READ (1,*) NIPT
      IF ( NIPT .GE. 1 ) THEN
      DO J = 1  , NIPT
      READ (1,*) I, XI(I), YI(I)
      END DO
      END IF
      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 )
      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*DY - RY*DX ) / (R*R) * W(IGAUSS)
      END DO
      GE =  C1 * DETJ * GE
      FE = -C1 * DETJ /DS * FE
      RETURN
      END
C
C
      SUBROUTINE INTE2 (INTEPT, ND, XP, YP, C1, XE, YE, SAI,W,
     *                   GEX, GEY, FEX, FEY )
      IMPLICIT REAL*8 ( A-H , O-Z )
      DIMENSION XE(ND), YE(ND), SAI(INTEPT), W(INTEPT)
      GEX = 0.D0
      GEY = 0.D0
      FEX = 0.D0
      FEY = 0.D0
      DX = XE(2) - XE(1)
      DY = YE(2) - YE(1)
      DS = DSQRT ( DX*DX + DY*DY )
      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)
      GEX = GEX + RX / (R*R) * W(IGAUSS)
      GEY = GEY + RY / (R*R) * W(IGAUSS)
C-------- INTEGRATION OF F(R)
      FEX = FEX + (-R*DY +2.D0*(RX*DY-RY*DX)*RX/R )/(R*R*R)*W(IGAUSS)
      FEY = FEY + ( R*DX +2.D0*(RX*DY-RY*DX)*RY/R )/(R*R*R)*W(IGAUSS)
      END DO
      GEX = -C1 * DETJ * GEX
      GEY = -C1 * DETJ * GEY
      FEX = -C1 * DETJ /DS * FEX
      FEY = -C1 * DETJ /DS * FEY
      RETURN
      END
C
C
      SUBROUTINE DOMAIN ( INTEPT,ND,MXE,MXI,C1,NIPT,NE,
     *                   SAI,W,XI,YI,NODEX,X,Y,H,QN,HI, CI, XE, YE )
      IMPLICIT REAL*8 ( A-H , O-Z )
      DIMENSION NODEX(MXE,ND), SAI(INTEPT), W(INTEPT),
     *  XI(MXI),YI(MXI), X(MXE),Y(MXE),H(MXE),QN(MXE),
     *  XE(ND),YE(ND),HI(MXI), CI(MXI)
C
      IF ( NIPT .EQ. 0 ) RETURN
C
      DO INSIDE = 1 , NIPT
      XP = XI(INSIDE)
      YP = YI(INSIDE)
      SUM = 0.D0
      C = 0.D0
      DO IEL = 1 , NE
      DO I = 1 , ND
      XE(I) = X( NODEX(IEL,I) )
      YE(I) = Y( NODEX(IEL,I) )
      END DO
      CALL INTE (INTEPT,ND,XP,YP,C1,XE,YE,SAI,W, G,F )
      C = C + F
      SUM = SUM + F*H(IEL) - G*QN(IEL) 
      END DO
      CI(INSIDE) = C
      HI(INSIDE) = SUM
      END DO
      RETURN
      END
C
C
      SUBROUTINE DOMAIN2 ( INTEPT,ND,MXE,MXI,C1,NIPT,NE,
     *                   SAI,W,XI,YI,NODEX,X,Y,H,QN,QX,QY,XE,YE )
      IMPLICIT REAL*8 ( A-H , O-Z )
      DIMENSION NODEX(MXE,ND), SAI(INTEPT), W(INTEPT),
     *  XI(MXI),YI(MXI), X(MXE),Y(MXE),H(MXE),QN(MXE),
     *  XE(ND),YE(ND),QX(MXI), QY(MXI)
C
      IF ( NIPT .EQ. 0 ) RETURN
C
      DO INSIDE = 1 , NIPT
      XP = XI(INSIDE)
      YP = YI(INSIDE)
      SUMX = 0.D0
      SUMY = 0.D0
      DO IEL = 1 , NE
      DO I = 1 , ND
      XE(I) = X( NODEX(IEL,I) )
      YE(I) = Y( NODEX(IEL,I) )
      END DO
      CALL INTE2 (INTEPT,ND,XP,YP,C1,XE,YE,SAI,W, GX,GY,FX,FY )
      SUMX = SUMX - FX*H(IEL) + GX*QN(IEL) 
      SUMY = SUMY - FY*H(IEL) + GY*QN(IEL) 
      END DO
      QX(INSIDE) = SUMX
      QY(INSIDE) = SUMY
      END DO
      RETURN
      END
C
C
      SUBROUTINE ECHOSOLN( MXE,MXI,NE,ND,NIPT,NODEX,X,Y,IELTYPE,H,QN,
     *                  BV, XI, YI, HI,CI, C, QX, QY )
      IMPLICIT REAL*8 ( A-H , O-Z )
      DIMENSION X(MXE), Y(MXE), XI(MXI), YI(MXI), HI(MXI),CI(MXI), 
     *  NODEX(MXE,ND), IELTYPE(MXE), H(MXE), QN(MXE), C(MXE),
     *  QX(MXI), QY(MXI), BV(MXE)
      CHARACTER*9 BC
      OPEN ( unit=2, FILE='SOLUTION.BEM', STATUS='UNKNOWN' )
C---------INPUT COORDINATES AND B.C.
      WRITE (2,*) "ELEMENT# X1 Y1 X2 Y2 BC BV"
      DO IEL = 1 , NE
      X1 = X(NODEX(IEL,1))
      Y1 = Y(NODEX(IEL,1))
      X2 = X(NODEX(IEL,2))
      Y2 = Y(NODEX(IEL,2))
      IF ( IELTYPE(IEL) .EQ. 1 ) THEN
      BC ='DIRICHLET'
      ELSE
      BC ='NUEMANN'
      END IF
      WRITE (2,*) IEL, X1, Y1, X2, Y2, BC, BV(IEL)
      END DO
C--------- SOLUTIONS ON BOUNADRY
      DIST = 0.D0
      WRITE (2,*)
      WRITE (2,*)"ELEMENT# X-CENTER Y-CENTER DIST POTENTIAL DPDN C"
      DO IEL = 1 , NE
      X1 = X(NODEX(IEL,1))
      Y1 = Y(NODEX(IEL,1))
      X2 = X(NODEX(IEL,2))
      Y2 = Y(NODEX(IEL,2))
      XCENTER = (X1+X2)/2.D0
      YCENTER = (Y1+Y2)/2.D0
      DL = DSQRT((X2-X1)**2 + (Y2-Y1)**2)
      DIST = DIST + DL/2.D0
      WRITE(2,*) IEL,XCENTER,YCENTER,DIST,H(IEL),QN(IEL),C(IEL)
      DIST = DIST + DL/2.D0
      END DO
C--------- INTERNAL VALUES
      WRITE (2,*)
      WRITE (2,*)"INTERNAL-POINT# XI YI POTENTIAL DPDX DPDY CI"
      DO I = 1 , NIPT
      WRITE(2,*) I, XI(I), YI(I), HI(I), QX(I), QY(I), CI(I)
      END DO
      CLOSE (2)
      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. 6 ) STOP'N>6'
      GO TO ( 99, 20, 30, 40, 50, 60 ) , N
   99 STOP
   20 SAI(1) = DSQRT(3.D0)/3.D0
      W(1) = 1.D0
      GO TO 88
   30 SAI(1) = DSQRT(15.D0)/5.D0
      SAI(2) = 0.D0
      W(1) = 5.D0/ 9.D0
      W(2) = 8.D0/ 9.D0
      GO TO 88
   40 SAI(1) = 0.33998104358485D0
      SAI(2) = 0.86113631159405D0
        W(1) = 0.65214515486254D0
        W(2) = 0.34785484513745D0
      GO TO 88
   50 SAI(1) = 0.90617984593866D0
      SAI(2) = 0.53846931010568D0
      SAI(3) = 0.D0
        W(1) = 0.23692688505619D0
        W(2) = 0.47862867049937D0
        W(3) = 5.12D0 / 9.D0
      GO TO 88
   60 SAI(1) = 0.23861918608320D0
      SAI(2) = 0.66120938646626D0
      SAI(3) = 0.93246951420315D0
        W(1) = 0.46791393457269D0
        W(2) = 0.36076157304814D0
        W(3) = 0.17132449237917D0
   88 NN = N / 2
      DO 11 I = 1 , NN
      J = N - I + 1
      SAI(J) = - SAI(I)
      W(J) = W(I)
   11 CONTINUE
      RETURN
      END
C
C
      SUBROUTINE SYSTEM ( MXN , N , A , C )
      IMPLICIT REAL*8 ( A-H , O-Z )
      DIMENSION A (MXN , MXN ) , C (MXN)
      N1 = N - 1
      DO 40 K = 1, N1
      L = K + 1
      DO 20 I = L , N
      P = A (I,K) / A (K,K)
      DO 30 J = L , N
   30 A (I,J) = A (I,J) - P * A ( K , J )
      C ( I ) = C ( I) - P * C ( K )
   20 CONTINUE
   40 CONTINUE
C---- BACK SUBSTITUTION
      C (N) = C (N) / A (N,N)
      DO 60 K = 1 , N1
      I = N - K
      L = I + 1
      P = C ( I )
      DO 50 J = L , N
   50 P = P - C (J) * A (I,J)
      C ( I ) = P / A (I,I)
   60 CONTINUE
      RETURN
      END