PROGRAM BEM8LINQWARP
C======================================================================
C BOUNDARY ELEMENT METHOD
C APPLIED TO
C LAPLACE EQUATION ( WARPING EQUATION )
C ELEMENT TYPE: LINEAR ELEMENT
C PROGRAMMED BY EIJI FUKUMORI // SUNY AT BUFFALO // 1984 SPRING
C 2 DEC. 2024
C======================================================================
IMPLICIT REAL*8 ( A-H , O-Z )
INCLUDE 'PARAM.DAT'
CCCCCCCCCCCCC PARAMETER (MXE=60, MXN=MXE, MXI=20,INTEPT=6, ND=2)
DIMENSION XE(ND),YE(ND),GE(ND),FE(ND),
* SAI(INTEPT),W(INTEPT),
* A(MXN,MXN),
* NODEX(MXE,ND),IELTYPE(MXE),BV(MXE,ND),
* C(MXN),
* X(MXN),Y(MXN),
* QN(MXN),H(MXN),RHS(MXN),NDTYPE(MXN),
* XI(MXI),YI(MXI), HI(MXI), CI(MXI)
C========================= INITIAL SET-UPS ============================
C
PI = 4.D0*DATAN( 1.D0)
C1 = - 1.D0/ ( 2.D0*PI )
CALL GRULE ( INTEPT, SAI, W )
C========================= READING IN DATA ============================
C
CALL INPUT (ND,MXE,MXI,NE,NODEX,IELTYPE,BV,X,Y,NIPT,XI,YI,
* GM, THETA )
C ...TYPE(I)=1 ---> TEMPERATURE PRESCRIBED
C ...TYPE(I)=2 ---> HEAT FLUX PRESCRIBED
NNODE = NE
C======================================================================
C================= ANALYSIS ON NODAL CONDITION ========================
C
CALL NDANALYS (ND,MXE,MXN,NE,IELTYPE,NODEX,NDTYPE )
C====================== DIRICHLET B.C. TO H(I) ========================
C
CALL BVARRANG ( ND, MXE,MXN, NE,NODEX,IELTYPE, BV, H )
C=============== FORMATION OF MATRIX A AND VECTOR RHS ==============
C
CALL MTXFORM (ND,MXE,MXN,C1,INTEPT,SAI,W,NE,NODEX,X,Y,
* NDTYPE, XE,YE, GE,FE, C ,IELTYPE,BV,RHS,H, A )
C===================== READY TO SOLVE A . X = C ===================
C
CALL SYSTEM ( MXE , NNODE, A , RHS )
C====================== SORTING SOLUTION ==============================
C
CALL SORTSOLN (ND,MXE,MXN,NE,IELTYPE,NODEX,RHS,NDTYPE,
* H,QN,BV )
C================== INTERNAL POINTS =================================
C
CALL DOMAIN ( INTEPT,ND,MXE,MXN,MXI,C1,NIPT,NE,
* GE,FE, SAI,W,XI,YI,NODEX,X,Y,H,BV,HI, CI,XE,YE )
C================= PRINTING RESULTS =================================
CALL ECHOSOLN( MXE,MXN,MXI,NE,ND,NIPT,NODEX,X,Y,IELTYPE,
* H, BV, XI, YI, HI,CI, C, THETA )
STOP 'NORMAL TERMINATION'
END
C
C
SUBROUTINE NDANALYS (ND,MXE,MXN,NE,IELTYPE,NODEX,NDTYPE )
DIMENSION IELTYPE(MXE), NDTYPE(MXN), NODEX(MXE,ND)
C------- THIS EVALUATES NODAL BOUNDARY CONDITION -------
C NDTYPE(NODE)=1 IMPLIES H IS KNOWN AT THE NODE.
C NDTYPE(NODE)=2 IMPLIES QN IS KNOWN AT THE NODE.
NNODE = NE
DO I = 1 , NNODE
NDTYPE(I) = 0
END DO
C------- CORNER(NODAL) BOUNDARY CONDITION EVALUATION -------
DO IEL = 1 , NE
DO J = 1 , ND
NDTYPE(NODEX(IEL,J)) = NDTYPE(NODEX(IEL,J)) + IELTYPE(IEL)
END DO
END DO
C-------
DO I = 1 , NNODE
IF ( NDTYPE(I) .EQ. 2 ) J = 1
IF ( NDTYPE(I) .EQ. 3 ) J = 3
IF ( NDTYPE(I) .EQ. 4 ) J = 2
NDTYPE(I) = J
END DO
RETURN
END
C
C
SUBROUTINE SORTSOLN (ND,MXE,MXN,NE,IELTYPE,NODEX,RHS,NDTYPE,
* H,QN,BV )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION IELTYPE(MXE),RHS(MXN),QN(MXN),H(MXN),NDTYPE(MXN),
* BV(MXE,ND),NODEX(MXE,ND)
NNODE = NE
DO I = 1 , NNODE
IF ( NDTYPE(I) .EQ. 2 ) THEN
H (I) = RHS(I)
ELSE
QN(I) = RHS(I)
END IF
END DO
C------- REDISTRIBUTION OF QN(I) INTO BV(I,J)
DO IEL = 1 , NE
IF ( IELTYPE(IEL) .EQ. 1 ) THEN
DO J = 1 , ND
NODE = NODEX(IEL,J)
BV(IEL,J) = QN(NODE)
END DO
END IF
END DO
RETURN
END
C
C
SUBROUTINE BVARRANG ( ND, MXE,MXN, NE,NODEX,IELTYPE, BV, H )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION IELTYPE(MXE),BV(MXE,ND),H(MXN),NODEX(MXE,ND)
C------ CONVERT TEMPRERATURE ON ELEMENT TO NODAL TEMPERATURE -
DO I = 1 , NE
DO J = 1 , ND
IF ( IELTYPE(I) .EQ. 1 ) THEN
H(NODEX(I,J)) = BV(I,J)
END IF
END DO
END DO
RETURN
END
C
C
SUBROUTINE MTXFORM (ND,MXE,MXN,C1,INTEPT,SAI,W,NE,NODEX,X,Y,
* NDTYPE, XE,YE, GE,FE, C ,IELTYPE,BV,RHS,H, A )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION C(MXN),X(MXN),Y(MXN),GE(ND),FE(ND),
* XE(ND),YE(ND),SAI(INTEPT),W(INTEPT) , NODEX(MXE,ND),
* IELTYPE(MXE),BV(MXE,ND), H(MXN),RHS(MXN),A(MXN,MXN),
* NDTYPE(MXN)
C=============== FORMATION OF MATRIX G,F AND VECTOR C
C---------- CLEAR MATRIX A(I,J) AND RHS(I) AND C(I)
NNODE = NE
DO I = 1 , NNODE
C(I) = 0.D0
RHS(I) = 0.D0
DO J = 1 , NNODE
A(I,J) = 0.D0
END DO
END DO
C ....... (XP,YP) = COORDINATES OF OBSERVATION POINT.
DO ISOURCE = 1 , NNODE
XP = X(ISOURCE)
YP = Y(ISOURCE)
DO IEL = 1 , NE
DO I = 1 , ND
XE(I) = X(NODEX(IEL,I))
YE(I) = Y(NODEX(IEL,I))
END DO
C------------------ INTEGRAL ON AN ELEMENT
NDDIFF1 = ISOURCE-NODEX(IEL,1)
NDDIFF2 = ISOURCE-NODEX(IEL,2)
IF ( NDDIFF1*NDDIFF2 .EQ. 0 ) THEN
CALL FINE ( ND,NDDIFF2, 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
DO L = 1 , ND
ICURREN = NODEX(IEL,L)
IF (IELTYPE(IEL) .EQ. 1 ) THEN
RHS(ISOURCE) = RHS(ISOURCE) + FE(L) * H(ICURREN)
A(ISOURCE,ICURREN) = A(ISOURCE,ICURREN) + GE(L)
END IF
IF (IELTYPE(IEL) .EQ. 2 ) THEN
RHS(ISOURCE) = RHS(ISOURCE) - GE(L) * BV(IEL,L)
IF ( NDTYPE(ICURREN) .EQ. 2 ) THEN
A(ISOURCE,ICURREN) = A(ISOURCE,ICURREN) - FE(L)
ELSE
RHS(ISOURCE) = RHS(ISOURCE) + FE(L) * H(ICURREN)
END IF
END IF
C(ISOURCE) = C(ISOURCE) + FE(L)
END DO
END DO
C.............FREE TERM
IF (NDTYPE(ISOURCE) .EQ. 2 ) THEN
A(ISOURCE,ISOURCE) = A(ISOURCE,ISOURCE) + C(ISOURCE)
ELSE
RHS(ISOURCE) = RHS(ISOURCE) - C(ISOURCE) * H(ISOURCE)
END IF
END DO
RETURN
END
C
C
SUBROUTINE INPUT (ND,MXE,MXI,NE,NODEX,IELTYPE,BV,X,Y,
* NIPT,XI,YI, GM, THETA )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION X(MXE),Y(MXE),IELTYPE(MXE),BV(MXE,ND),NODEX(MXE,ND),
* XI(MXI),YI(MXI)
OPEN ( 1, FILE='BEM1.DAT', STATUS='OLD' )
READ (1,*) GM, THETA
READ (1,*) NE
DO IEL = 1 , NE
READ (1,*) I,(NODEX(I,J),J=1,ND),IELTYPE(I), (BV(I,J),J=1,ND)
END DO
NNODE = NE
DO I = 1 , NNODE
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,NDDIFF2, C1, XE, YE, GE, FE )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION XE(ND),YE(ND), GE(ND),FE(ND)
FE(1) = 0.D0
FE(2) = 0.D0
DX = XE(2) - XE(1)
DY = YE(2) - YE(1)
DS = DSQRT ( DX*DX + DY*DY )
GE(1) = C1 * (DS/2.D0) *( DLOG(DS) - 1.5D0 )
GE(2) = C1 * (DS/2.D0) *( DLOG(DS) - 0.5D0 )
IF ( NDDIFF2 .EQ. 0 ) THEN
TEMP = GE(1)
GE(1) = GE(2)
GE(2) = TEMP
END IF
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(ND),FE(ND)
DO I = 1 , ND
GE(I) = 0.D0
FE(I) = 0.D0
END DO
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--------GAUSS INTEGRATION
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 )
SF1 = 0.5D0 * ( 1.D0 - SAI(IGAUSS) )
SF2 = 0.5D0 * ( 1.D0 + SAI(IGAUSS) )
C-------- INTEGRATION OF G(R)
TEMP = DLOG(R) * W(IGAUSS)
GE(1) = GE(1) + TEMP * SF1
GE(2) = GE(2) + TEMP * SF2
C-------- INTEGRATION OF F(R)
TEMP = (RX*DY-RY*DX) / (R*R) * W(IGAUSS)
FE(1) = FE(1) + TEMP * SF1
FE(2) = FE(2) + TEMP * SF2
END DO
GE(1) = C1 * DETJ * GE(1)
GE(2) = C1 * DETJ * GE(2)
FE(1) = -C1 * DETJ /DS * FE(1)
FE(2) = -C1 * DETJ /DS * FE(2)
RETURN
END
C
C
SUBROUTINE DOMAIN ( INTEPT,ND,MXE,MXN,MXI,C1,NIPT,NE,
* GE,FE, SAI,W,XI,YI,NODEX,X,Y,H,BV,HI, CI,XE,YE )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION NODEX(MXE,ND), SAI(INTEPT), W(INTEPT),
* XI(MXI),YI(MXI), X(MXN),Y(MXN),H(MXN), BV(MXE,ND),
* XE(ND),YE(ND),HI(MXI),CI(MXI),
* GE(ND),FE(ND)
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, GE,FE )
DO J = 1 , ND
C = C + FE(J)
SUM = SUM + FE(J)*H(NODEX(IEL,J)) - GE(J)*BV(IEL,J)
END DO
END DO
CI(INSIDE) = C
HI(INSIDE) = SUM
END DO
RETURN
END
C
C
SUBROUTINE ECHOSOLN( MXE,MXN,MXI,NE,ND,NIPT,NODEX,X,Y,IELTYPE,
* H, BV, XI, YI, HI,CI, C , THETA )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION X(MXN), Y(MXN), XI(MXI), YI(MXI), HI(MXI),CI(MXI),
* NODEX(MXE,ND), IELTYPE(MXE), H(MXN), BV(MXE,ND), C(MXN)
CHARACTER*9 BC
OPEN ( 1, FILE='SOLUTION.BEM', STATUS='UNKNOWN' )
C---------INPUT COORDINATES AND B.C.
WRITE (1,*) '========== STEADY HEAT TRANSFER ANALYSIS =========='
WRITE (1,*) '=========== BY BEM USING LINEAR ELEMENT ==========='
WRITE (1,*) '== ALONG WITH STRAIGHT FORMATION OF [A]{X}={RHS} =='
WRITE (1,*)
WRITE (1,*) ' ==== ELEMENT INFORMATION ===='
WRITE (1,*) " ELEMENT# NODE(1) NODE(2) BC B-VALUE(1) B-VALUE(2)"
DO IEL = 1 , NE
IF ( IELTYPE(IEL) .EQ. 1 ) THEN
BV1 = H(NODEX(IEL,1))
BV2 = H(NODEX(IEL,2))
BC ='DIRICHLET'
END IF
IF ( IELTYPE(IEL) .EQ. 2 ) THEN
BV1 = BV(IEL,1)
BV2 = BV(IEL,2)
BC ='NUEMANN'
END IF
WRITE (1,*) IEL, NODEX(IEL,1), NODEX(IEL,2), BC, BV1, BV2
END DO
C
WRITE(1,*)
WRITE(1,*) ' == NODAL INFOMATION AND SOLUTION OF TEMPERATURE =='
WRITE(1,*) "NODE# X-COORD Y-COORD S-COORD PSI W(X,Y) FREE-TERM"
NNODE = NE
S = 0.D0
I = 1
WRITE (1,*) I, X(I) , Y(I) , S, H(I) , THETA*H(I), C(I)
DO I = 2 , NNODE
S = S + DSQRT ( (X(I)-X(I-1))**2 + (Y(I)-Y(I-1))**2 )
WRITE (1,*) I, X(I) , Y(I) , S, H(I) , THETA*H(I), C(I)
END DO
I = 1
S = S + DSQRT ( (X(I)-X(NNODE))**2 + (Y(I)-Y(NNODE))**2 )
WRITE (1,*) I, X(I) , Y(I) , S, H(I) , THETA*H(I), C(I)
C
WRITE(1,*)
WRITE(1,*) '=== SOLUTION OF HEAT FLUX ==='
WRITE(1,*)" ELEMENT# HEAT-FLUX(1) HEAT-FLUX(2)"
DO IEL = 1 , NE
WRITE(1,*) IEL, BV(IEL,1), BV(IEL,2)
END DO
C
WRITE(1,*)
WRITE(1,*) ' === TEMPERATURE AT INTERNAL POINTS ==='
WRITE (1,*)"POINT# X-COORD Y-COORD TEMPERATURE FREE-TERM"
DO I = 1 , NIPT
WRITE(1,*) I, XI(I), YI(I), HI(I), CI(I)
END DO
CLOSE (1)
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 K = 1, N1
L = K + 1
DO I = L , N
P = A (I,K) / A (K,K)
IF ( P .NE. 0. ) THEN
DO J = L , N
A (I,J) = A (I,J) - P * A ( K , J )
END DO
C ( I ) = C ( I) - P * C ( K )
END IF
END DO
END DO
C---- BACK SUBSTITUTION
C (N) = C (N) / A (N,N)
DO K = 1 , N1
I = N - K
L = I + 1
P = C ( I )
DO J = L , N
P = P - C (J) * A (I,J)
END DO
C ( I ) = P / A (I,I)
END DO
RETURN
END
C
C
SUBROUTINE GRULE ( N , SAI , W )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION SAI(N) , W(N), TABLESAI(8,4), TABLEW(8,4)
IF ( N .LT. 1 ) STOP'N<1'
IF ( N .GT. 8 ) STOP'N>8'
C-------- 1-POINT GAUSS-LEGENDRE
TABLESAI(1,1) = 0.D0
TABLEW(1,1) = 2.D0
C-------- 2-POINT GAUSS-LEGENDRE
TABLESAI(2,1) = DSQRT(3.D0)/3.D0
TABLEW(2,1) = 1.D0
C-------- 3-POINT GAUSS-LEGENDRE
TABLESAI(3,1) = DSQRT(15.D0)/5.D0
TABLESAI(3,2) = 0.D0
TABLEW(3,1) = 5.D0/ 9.D0
TABLEW(3,2) = 8.D0/ 9.D0
C-------- 4-POINT GAUSS-LEGENDRE
TABLESAI(4,1) = 0.33998104358485626480D0
TABLESAI(4,2) = 0.86113631159405257522D0
TABLEW(4,1) = 0.6521451548625461426D0
TABLEW(4,2) = 0.34785484513745385737D0
C-------- 5-POINT GAUSS-LEGENDRE
TABLESAI(5,1) = 0.90617984593866399279D0
TABLESAI(5,2) = 0.53846931010568309103D0
TABLESAI(5,3) = 0.D0
TABLEW(5,1) = 0.23692688505618908751D0
TABLEW(5,2) = 0.47862867049936646804D0
TABLEW(5,3) = 5.12D0 / 9.D0
C-------- 6-POINT GAUSS-LEGENDRE
TABLESAI(6,1) = 0.23861918608319690863D0
TABLESAI(6,2) = 0.66120938646626451366D0
TABLESAI(6,3) = 0.93246951420315202781D0
TABLEW(6,1) = 0.46791393457269104738D0
TABLEW(6,2) = 0.36076157304813860756D0
TABLEW(6,3) = 0.17132449237917034504D0
C-------- 7-POINT GAUSS-LEGENDRE
TABLESAI(7,1) = 0.94910791234275852452D0
TABLESAI(7,2) = 0.74153118559939443986D0
TABLESAI(7,3) = 0.40584515137739716690D0
TABLESAI(7,4) = 0.D0
TABLEW(7,1) = 0.12948496616886969327D0
TABLEW(7,2) = 0.27970539148927666790D0
TABLEW(7,3) = 0.38183005050511894495D0
TABLEW(7,4) = 0.41795918367346938775D0
C-------- 8-POINT GAUSS-LEGENDRE
TABLESAI(8,1) = 0.96028985649753623168D0
TABLESAI(8,2) = 0.79666647741362673959D0
TABLESAI(8,3) = 0.52553240991632898581D0
TABLESAI(8,4) = 0.18343464249564980493D0
TABLEW(8,1) = 0.10122853629037625915D0
TABLEW(8,2) = 0.22238103445337447054D0
TABLEW(8,3) = 0.31370664587788728733D0
TABLEW(8,4) = 0.36268378337836198296D0
C----- SAI BETWEEN 0 AND +1
NFRONT = N / 2
NREAR = N - NFRONT
DO I = 1, NREAR
SAI(I) = TABLESAI(N,I)
W(I) = TABLEW(N,I)
END DO
C----- FOR A CASE OF N=1
IF ( NFRONT .EQ. 0 ) THEN
RETURN
END IF
C----- SAI BETWEEN -1 AND 0
DO I = 1 , NFRONT
J = N - I + 1
SAI(J) = - TABLESAI(N,I)
W(J) = TABLEW(N,I)
END DO
RETURN
END