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