PROGRAM POSTPROCESS
C======================================================================
C POST PROCESS PROGRAM
C ------------- UNBOUNDED VERSION -------------
C BOUNDARY ELEMENT METHOD
C APPLIED TO
C LAPLACE EQUATION (STEADY HEAT EQUATION )
C ELEMENT TYPE: LINEAR ELEMENT
C PROGRAMMED BY EIJI FUKUMORI // NAGOYA // 2021DEC
C======================================================================
IMPLICIT REAL*8 ( A-H , O-Z )
PARAMETER (MXE=1000, MXN=MXE, MXI=1000,INTEPT=6, ND=2)
PARAMETER ( NDVEC=4,MXEVEC=5000000,MXNVEC=5500000 )
DIMENSION NODEX(MXE,ND),X(MXN),Y(MXN),H(MXN),
* XI(MXI),YI(MXI), HI(MXI), ANGLE(MXN), BV(MXE,ND)
DIMENSION NODEXVEC(MXEVEC,NDVEC), XCOORDVEC(MXNVEC)
DIMENSION YCOORDVEC(MXNVEC), HVEC(MXNVEC)
DIMENSION XE(ND),YE(ND),GE(ND),FE(ND),SAI(INTEPT),W(INTEPT)
C23456789-123456789-123456789-123456789-123456789-123456789-123456789-12
C========================= READING IN DATA ============================
C
CALL INPUT (ND,MXE,MXI,NE,NODEX,X,Y,NIPT,XI,YI,
* MXEVEC,NDVEC,MXNVEC,NEVEC,NNODEVEC,
* NODEXVEC,XCOORDVEC,YCOORDVEC )
NNODE = NE
NC = NNODE/2
PI = 4.D0*DATAN(1.D0)
CALL READSOL ( ND,MXE,MXN,MXI,NNODE,NIPT, H, HI, BV )
OPEN ( 3, FILE='BOUNDARY.SOL', STATUS='UNKNOWN' )
WRITE (3,*) 'ANGLE[RADIAN] POTENTIAL FLUX'
YMAX = Y(1)
YMIN = Y(1)
DO I = 2, NC
IF ( Y(I) .GT. YMAX ) YMAX = Y(I)
IF ( Y(I) .LT. YMIN ) YMIN = Y(I)
END DO
YC1 = (YMAX + YMIN) /2.D0
CCCCCCCCC WRITE (3,*) 'YC1=', YC1
DO I = 1 , NC
YNODE = Y(I) - YC1
ANGLE (I) = ANGLERD (PI,X(I),YNODE)
WRITE (3,*) ANGLE (I), H(I), BV(I,1)
END DO
WRITE (3,*)
BTINTE = 0.D0
DO IEL = 1 , NE/2
I = NODEX(IEL,1)
J = NODEX(IEL,2)
DX = X(J) - X(I)
DY = Y(J) - Y(I)
DS = DSQRT (DX*DX+DY*DY)
BTINTE = BTINTE -DS*(BV(IEL,1)+BV(IEL,2))/2.D0
END DO
C---------------------------------------------------------------
YMAX = Y(1+NC)
YMIN = Y(1+NC)
DO I = NC+2, NNODE
IF ( Y(I) .GT. YMAX ) YMAX = Y(I)
IF ( Y(I) .LT. YMIN ) YMIN = Y(I)
END DO
YC2 = (YMAX + YMIN) /2.D0
DO I = NC+1 , NNODE
YNODE = Y(I) - YC2
ANGLE (I) = ANGLERD (PI,X(I),YNODE)
WRITE (3,*) ANGLE (I), H(I), BV(I,1)
END DO
CLOSE (3)
OPEN ( 4, FILE='DELTAP.SOL', STATUS='UNKNOWN' )
DELTAP = DABS (H(1) - H(NC+1))
DO I = 2, NC
DELTAP = DELTAP + DABS (H(I) - H(I+NC))
END DO
DELTAP = DELTAP/NC
WRITE (4,*) 'DELTA-P=',DELTAP
WRITE (4,*) 'LINE INTEGRAL OF AMPERE OR GAUSS(BTINTE) =', BTINTE
WRITE (4,*) '1/(BTINTE/DELTA-P)=',1.D0/(BTINTE/DELTAP)
CLOSE (4)
C23456789-123456789-123456789-123456789-123456789-123456789-123456789-12
C============================================================================
C
C VECTOR PLOT AND CONTOUR PLOT ROUTINE
C
C============================================================================
C1 = - 1.D0/ ( 8.D0* DATAN( 1.D0) )
CALL GRULE ( INTEPT, SAI, W )
CALL PLOT ( INTEPT,ND,MXE,MXN,C1,NE,GE,FE,SAI,W,XE,YE,
* NODEX,X,Y,H,BV, MXEVEC,NDVEC,MXNVEC,NEVEC,NNODEVEC,
* XCOORDVEC,YCOORDVEC, HVEC )
CALL CONTOUR (MXEVEC,MXNVEC,NDVEC,NEVEC,NNODEVEC,NODEXVEC,
* XCOORDVEC,YCOORDVEC, HVEC )
CALL VECTOR ( NDVEC,MXEVEC,MXNVEC,NEVEC,NNODEVEC,NODEXVEC,
* XCOORDVEC,YCOORDVEC, HVEC )
STOP 'NORMAL TERMINATION'
END
C
C
C23456789-123456789-123456789-123456789-123456789-123456789-123456789-12
SUBROUTINE PLOT ( INTEPT,ND,MXE,MXN,C1,NE,GE,FE,SAI,W,XE,YE,
* NODEX,X,Y,H,BV, MXEVEC,NDVEC,MXNVEC,NEVEC,NNODEVEC,
* XCOORDVEC,YCOORDVEC, HVEC )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION NODEX(MXE,ND), SAI(INTEPT), W(INTEPT), X(MXN),Y(MXN),
* H(MXN), BV(MXE,ND),XE(ND),YE(ND), GE(ND),FE(ND)
DIMENSION XCOORDVEC(MXNVEC), YCOORDVEC(MXNVEC), HVEC(MXNVEC)
C
C
DO INSIDE = 1 , NNODEVEC
XP = XCOORDVEC(INSIDE)
YP = YCOORDVEC(INSIDE)
SUM = 0.D0
C = 1.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
IC = ABS(C)+0.01
IF ( IC .NE. 1 ) THEN
IF ( YP .GT. 0.D0 ) SUM = H(NODEX(1,1))
IF ( YP .LT. 0.D0 ) SUM = H(NODEX(NE,1))
END IF
HVEC(INSIDE) = SUM
END DO
RETURN
END
C
C
SUBROUTINE VECTOR ( NDVEC,MXEVEC,MXNVEC,NEVEC,NNODEVEC,NODEXVEC,
* XCOORDVEC,YCOORDVEC, HVEC )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION NODEXVEC(MXEVEC,NDVEC),XCOORDVEC(MXNVEC),
* YCOORDVEC(MXNVEC),HVEC(MXNVEC)
C-------XCOORD AND YCOORD ARE COORDINATES OF CONDUCTOR------
C-------XCOORDVEC AND YCOORDVEC ARE COORDINATES OF observation point------
C=================================================================
C23456789-123456789-123456789-123456789-123456789-123456789-123456789-12
C----------------- COMPUTATION OF B VECTOR--------------------
IUNIT = 1
OPEN ( IUNIT, FILE='VECTORB.SOL', STATUS='UNKNOWN' )
JUNIT = 2
OPEN ( JUNIT, FILE='VECTORE.SOL', STATUS='UNKNOWN' )
ALPHA = 3.2D0
DO IEL = 1 , NEVEC
I = NODEXVEC(IEL,1)
J = NODEXVEC(IEL,2)
K = NODEXVEC(IEL,3)
L = NODEXVEC(IEL,4)
XP = (XCOORDVEC(I)+XCOORDVEC(J)+XCOORDVEC(K)+XCOORDVEC(L))/4.D0
YP = (YCOORDVEC(I)+YCOORDVEC(J)+YCOORDVEC(K)+YCOORDVEC(L))/4.D0
DX = XCOORDVEC(J) - XCOORDVEC(I)
DY = YCOORDVEC(K) - YCOORDVEC(J)
DADX = ( (HVEC(J)-HVEC(I))/DX+(HVEC(K)-HVEC(L))/DX )/2.D0
DADY = ( (HVEC(K)-HVEC(J))/DY+(HVEC(L)-HVEC(I))/DY )/2.D0
BX = DADY
BY = - DADX
CALL VECPLT ( IUNIT, XP, YP, BX, BY, ALPHA )
EX = - DADX
EY = - DADY
CALL VECPLT ( JUNIT, XP, YP, EX, EY, ALPHA )
END DO
CLOSE (IUNIT)
CLOSE (JUNIT)
RETURN
END
C
C
SUBROUTINE VECPLT ( IUNIT, X0, Y0, U, V, FACT )
IMPLICIT REAL*8 ( A-H , O-Z )
DATA AL, BETA / 0.65D0, 0.15D0 /
DX = U * FACT
DY = V * FACT
WRITE(IUNIT,*) X0, Y0
WRITE(IUNIT,*) X0+ AL*DX , Y0+ AL*DY
RNX = DY
RNY = -DX
X = X0 + AL*DX + BETA*RNX
Y = Y0 + AL*DY + BETA*RNY
WRITE (IUNIT,*) X, Y
WRITE(IUNIT,*) X0+DX , Y0+DY
X = X0 + AL*DX - BETA*RNX
Y = Y0 + AL*DY - BETA*RNY
WRITE (IUNIT,*) X, Y
WRITE(IUNIT,*) X0+ AL*DX , Y0+ AL*DY
WRITE(IUNIT,*)
RETURN
END
C
C
SUBROUTINE CONTOUR (MXEVEC,MXNVEC,NDVEC,NEVEC,NNODEVEC,NODEXVEC,
* XCOORDVEC,YCOORDVEC, P )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION NODEXVEC(MXEVEC,NDVEC), XCOORDVEC(MXNVEC),
* YCOORDVEC(MXNVEC), P(MXNVEC),
* S(4), B(2,4)
OPEN ( 1, FILE='CONTOUR.SOL', STATUS='UNKNOWN' )
CUTLIMIT = 1.0D0
NSTEP = 20
NSTEP = NSTEP/2*2 + 1
PPMAX = P(1)
PPMIN = P(1)
DO I = 2 , NNODEVEC
IF ( P(I) .GT. PPMAX ) PPMAX = P(I)
IF ( P(I) .LT. PPMIN ) PPMIN = P(I)
END DO
PPMAX = CUTLIMIT * PPMAX
PPMIN = CUTLIMIT * PPMIN
WRITE (*,*) PPMAX , PPMIN
DS = ( PPMAX - PPMIN ) / NSTEP
C
CCCCCCCCC IC = 0
DO IEL = 1 , NEVEC
DO I = 1 , NDVEC
B(1,I) = XCOORDVEC(NODEXVEC(IEL,I))
B(2,I) = YCOORDVEC(NODEXVEC(IEL,I))
S(I) = P(NODEXVEC(IEL,I))
END DO
IF ( (PPMIN .LT. 0.D0) .AND. (PPMAX.GT.0.D0) ) THEN
CALL PLTSAI ( DS, NSTEP, -NSTEP*DS, B, S )
CALL PLTSAI ( DS, 1, 0.D0, B, S )
CALL PLTSAI ( DS, NSTEP, DS, B, S )
ELSE
CALL PLTSAI ( DS, NSTEP, PPMIN, B, S )
END IF
END DO
CLOSE (1)
RETURN
END
C
C
SUBROUTINE INPUT (ND,MXE,MXI,NE,NODEX,X,Y,NIPT,XI,YI,
* MXEVEC,NDVEC,MXNVEC,NEVEC,NNODEVEC,
* NODEXVEC,XCOORDVEC,YCOORDVEC )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION X(MXE),Y(MXE),IELTYPE(MXE),NODEX(MXE,ND),
* XI(MXI),YI(MXI)
DIMENSION NODEXVEC(MXEVEC,NDVEC), XCOORDVEC(MXNVEC)
DIMENSION YCOORDVEC(MXNVEC)
OPEN ( 1, FILE='BEM1.DAT', STATUS='OLD' )
READ (1,*) FMU
READ (1,*) NE
DO IEL = 1 , NE
READ (1,*) I,(NODEX(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)
OPEN ( 5, FILE='VECTORCG.DAT', STATUS = 'UNKNOWN' )
READ(5,*) NEVEC
DO I = 1 , NEVEC
READ (5,*) IEL,(NODEXVEC(IEL,J),J=1,NDVEC)
END DO
READ(5,*) NNODEVEC
DO I = 1 , NNODEVEC
READ(5,*) NODE, XCOORDVEC(NODE), YCOORDVEC(NODE)
END DO
CLOSE (5)
RETURN
END
C
C
SUBROUTINE READSOL ( ND,MXE,MXN,MXI,NNODE,NIPT, H, HI, BV )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION HI(MXI), H(MXN), BV(MXE,ND)
C------------- READING A FILE FOR SOLUTION PROCESS --------------
OPEN ( 2, FILE='POSTPROC.DAT', STATUS='UNKNOWN' )
DO I = 1 , NNODE
READ (2,*) H(I)
END DO
DO I = 1 , NIPT
READ (2,*) HI(I)
END DO
NE = NNODE
DO I = 1 , NE
READ (2,*) IEL, BV(IEL,1), BV(IEL,2)
END DO
CLOSE (2)
RETURN
END
C
C
FUNCTION ANGLERD (PI,X,Y)
IMPLICIT REAL*8 ( A-H , O-Z )
IF ( Y .EQ. 0.D0 .AND. X .GT. 0.D0 ) THEN
ANGLERD =0.D0
RETURN
END IF
IF ( Y .GT. 0.D0 .AND. X .EQ. 0.D0 ) THEN
ANGLERD =0.5D0*PI
RETURN
END IF
IF ( Y .EQ. 0.D0 .AND. X .LT. 0.D0 ) THEN
ANGLERD =PI
RETURN
END IF
IF ( Y .LT. 0.D0 .AND. X .EQ. 0.D0 ) THEN
ANGLERD =1.5D0*PI
RETURN
END IF
IF ( X .NE. 0.D0 ) ANG = DATAN (DABS(Y)/DABS(X))
IF ( Y .GT. 0.D0 .AND. X .GT. 0.D0 ) THEN
ANGLERD =ANG
RETURN
END IF
IF ( Y .GT. 0.D0 .AND. X .LT. 0.D0 ) THEN
ANGLERD =PI-ANG
RETURN
END IF
IF ( Y .LT. 0.D0 .AND. X .LT. 0.D0 ) THEN
ANGLERD =PI+ANG
RETURN
END IF
IF ( Y .LT. 0.D0 .AND. X .GT. 0.D0 ) THEN
ANGLERD =2.D0*PI-ANG
RETURN
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 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 PLTSAI ( DS, NSTEP, START, CRD, SS )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION CRD(2,4), SS(4), X(4), Y(4), S(4)
IF ( NSTEP .EQ. 0 ) RETURN
X(3) = ( CRD(1,1) + CRD(1,2) + CRD(1,3) + CRD(1,4) ) / 4.
Y(3) = ( CRD(2,1) + CRD(2,2) + CRD(2,3) + CRD(2,4) ) / 4.
S(3) = ( SS(1) + SS(2) + SS(3) + SS(4) ) / 4.
SMAX = DMAX1 ( SS(1), SS(2), SS(3), SS(4) )
SMIN = DMIN1 ( SS(1), SS(2), SS(3), SS(4) )
DO 52 LEVEL = 1 , NSTEP
SXY = START + (LEVEL-1) * DS
IF ( (SMAX-SXY)*(SMIN-SXY) .LT. 0 ) THEN
DO 60 IEL = 1 , 4
X(1) = CRD(1,IEL)
Y(1) = CRD(2,IEL)
S(1) = SS(IEL)
IF ( IEL .EQ. 4 ) THEN
X(2) = CRD(1,1)
Y(2) = CRD(2,1)
S(2) = SS(1)
ELSE
X(2) = CRD(1,IEL+1)
Y(2) = CRD(2,IEL+1)
S(2) = SS (IEL+1)
ENDIF
X(4) = X(1)
Y(4) = Y(1)
S(4) = S(1)
K = 0
DO 70 ISG = 1 , 3
IF ( S(ISG ) .LT. SXY ) GO TO 30
IF ( S(ISG+1) .LT. SXY ) GO TO 40
GO TO 70
30 IF ( S(ISG+1) .LT. SXY ) GO TO 70
40 T = ( SXY - S(ISG) ) / ( S(ISG+1) - S(ISG) )
X0 = X(ISG+1)*T + (1.- T)*X(ISG)
Y0 = Y(ISG+1)*T + (1.- T)*Y(ISG)
IF ( K .EQ. 0 ) GO TO 71
CALL XDRAW ( X0, Y0 )
GO TO 60
71 CALL XMOVE ( X0, Y0 )
K = 1
70 CONTINUE
60 CONTINUE
ENDIF
52 CONTINUE
RETURN
END
C
C
SUBROUTINE XMOVE ( X0, Y0 )
IMPLICIT REAL*8 ( A-H , O-Z )
WRITE (1,*)
WRITE (1,*) X0, Y0
RETURN
END
C
C
SUBROUTINE XDRAW ( X0, Y0 )
IMPLICIT REAL*8 ( A-H , O-Z )
WRITE (1,*) X0, Y0
RETURN
END