PROGRAM WAVECOMP
C======================================================================
C COMPUTATION OF WAVE HEIGHT COEFFICIENT
C AT SELECTED POINTS IN INTERIOR OF COMPUTATIONAL DOMAIN
C
C NOTE: TOPOGRAPHY FUNCTION OF WAVE PROPERGATION MUST BE EVALUATED
C BY HLM8LINQ.FOR PRIOR TO THIS PROGRAM EXCECUTION
C
C ELEMENT TYPE FOR BOUNDARY: LINEAR ELEMENT
C ELEMENT TYPE FOR DOMAIN: FOUR-NODED ELEMENT
C INPUT FILES: BEM1.DAT & BOUNDARY.BIN
C PROGRAMMED BY EIJI FUKUMORI // 2025 JAN. 12
C======================================================================
IMPLICIT REAL*8 ( A-H , O-Z )
PARAMETER (MXE=150, MXN=MXE,INTEPT=6, ND=2 )
PARAMETER ( MXEI=50000, MXNI=50000,NDI=4 )
DIMENSION XE(ND),YE(ND),SAI(INTEPT),W(INTEPT),NODEX(MXE,ND),
* X(MXN),Y(MXN),WAVEHB(MXN)
DIMENSION XI(MXNI),YI(MXNI), NODEXI(MXEI,NDI),WAVEH(MXNI)
COMPLEX*16 BV(MXE,ND),GE(ND),FE(ND),H(MXN),C2
COMPLEX*16 HI(MXNI)
C========================= INITIAL SET-UPS ============================
C
PI = 4.D0 * DATAN( 1.D0)
C2 = 1.D0/ DCMPLX(0.,4.D0)
CALL GRULE ( INTEPT, SAI, W )
C========================= READING IN DATA ============================
C
CALL INPUT (OMEGA,ND,MXE,MXN,NE,NNODE,NODEX,BV,X,Y,H,WAVEHB)
C ...TYPE(I)=1 ---> TEMPERATURE PRESCRIBED
C ...TYPE(I)=2 ---> HEAT FLUX PRESCRIBED
C=============== EVALUATION OF ELEMENT SIZE AND OBJECT SIZES ==========
CALL SIZES (ND,MXE,MXN,NE,NNODE,NODEX,X,Y,DS,
* BODYX,BODYY,CENTERX, CENTERY )
C====== NODAL POINTS AND ELEMENTS CREATION FOR INTERIOR POINTS=========
CALL CREATE ( NDI,MXEI,MXNI,DS,BODYX,BODYY,CENTERX,CENTERY,
* NEI, NNODEI, NODEXI,XI,YI )
C================== INTERNAL POINTS =================================
C
CALL DOMAIN ( PI,OMEGA,INTEPT,ND,MXE,MXN,MXNI,C2,NNODEI,
* NE, GE,FE, SAI,W,XI,YI,NODEX,X,Y,H,BV,HI )
C============== COMPUTATION OF WAVE HEIGHT COEFFICIENT ================
CALL WAVEHT ( OMEGA, MXNI,NNODEI,HI, XI, WAVEH )
C========== MERGING BOUNDARY WAVE HEIGHT TO DOMAIN WAVE HEIGHT ========
CALL MERGE ( MXN,NNODE,X,Y,WAVEHB,MXNI,NNODEI,XI,YI,WAVEH )
C==================== CREATE FEM LIKE INPUT FILE ======================
CALL CRFILE ( NDI,MXEI,MXNI,NEI, NNODEI, NODEXI,XI,YI,WAVEH )
STOP 'NORMAL TERMINATION'
END
C
C
SUBROUTINE MERGE ( MXN,NNODE,X,Y,WAVEHB,MXNI,NNODEI,XI,YI,WAVEH )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION XI(MXNI),YI(MXNI),WAVEH(MXNI)
DIMENSION X(MXN),Y(MXN),WAVEHB(MXN)
DO I = 1 , NNODE
RMIN = DSQRT ( (XI(1)-X(I))**2 + (YI(1)-Y(I))**2 )
LOC = 1
DO J = 2 , NNODEI
R = DSQRT ( (XI(J)-X(I))**2 + (YI(J)-Y(I))**2 )
IF ( R .LT. RMIN ) THEN
LOC = J
RMIN = R
END IF
END DO
WAVEH(LOC) = WAVEHB(I)
XI(LOC) = X(I)
YI(LOC) = Y(I)
END DO
RETURN
END
C
C
SUBROUTINE CRFILE ( NDI,MXEI,MXNI,NEI,NNODEI,NODEXI,XI,YI,WAVEH )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION NODEXI(MXEI,NDI), XI(MXNI),YI(MXNI),WAVEH(MXNI)
OPEN (1,FILE='DOMAIN.BIN',STATUS='UNKNOWN',FORM='UNFORMATTED')
WRITE(1) NEI
WRITE(1) NNODEI
DO I = 1 , NEI
WRITE(1) (NODEXI(I,J),J=1,NDI)
END DO
WRITE(1) ( XI(I) , I=1 , NNODEI )
WRITE(1) ( YI(I) , I=1 , NNODEI )
WRITE(1) ( WAVEH(I) , I=1 , NNODEI )
CLOSE (1)
RETURN
END
C
C
SUBROUTINE CREATE ( NDI,MXEI,MXNI,DS,BODYX,BODYY,CENTERX,CENTERY,
* NEI, NNODEI, NODEXI,XI,YI )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION NODEXI(MXEI,NDI), XI(MXNI),YI(MXNI)
C------ IOFFXUP = SEE HOMEPAGE FOR EXPLANATION
C------ IOFFXDW = SEE HOMEPAGE FOR EXPLANATION
C------ IOFFY = SEE HOMEPAGE FOR EXPLANATION
DX = DS
DY = DS
IOFFXUP = 1
IOFFXDW = 1
IOFFY = 1
TLX = (IOFFXUP+1+IOFFXDW)*BODYX
TLY = (IOFFY +1+IOFFY )*BODYY
NEX = TLX/DS
NEY = TLY/DS
NEY = NEY/2*2
XORIGN = CENTERX - BODYX/2. - IOFFXUP*BODYX
YORIGN = CENTERY - (NEY/2)*DS
NDX = NEX + 1
NDY = NEY + 1
C=======================================================================
WRITE (*,*) 'DOMAIN LENGTH IN X (TLX) =',TLX
WRITE (*,*) 'DOMAIN LENGTH IN Y (TLY) =',TLY
WRITE (*,*) 'NUMBER OF ELEMENTS IN X (NEX) =', NEX
WRITE (*,*) 'NUMBER OF ELEMENTS IN Y (NEY) =', NEY
WRITE (*,*) 'INCREMENTAL LENGTH IN X (DX) =', DX
WRITE (*,*) 'INCREMENTAL LENGTH IN Y (DX) =', DY
WRITE (*,*) 'ORIGIN OF GLOBAL COORDINATES =',XORIGN,YORIGN
C=======================================================================
C ELEMENT CREATION
NEI = 0
DO J = 1 , NEX
DO I = 1 , NEY
NEI = NEI + 1
IF ( NEI .GT. MXEI ) STOP 'NEI > MXEI'
NODEXI(NEI,1) = NDY*(J-1) + I
NODEXI(NEI,2) = NODEXI(NEI,1) + NDY
NODEXI(NEI,3) = NODEXI(NEI,2) + 1
NODEXI(NEI,4) = NODEXI(NEI,1) + 1
END DO
END DO
C=======================================================================
C NODAL COORDINATE CREATION
NNODEI = 0
DO I = 1 , NDX
DO J = 1 , NDY
NNODEI = NNODEI + 1
IF ( NNODEI .GT. MXNI ) STOP 'NNODEI > MXNI'
XI(NNODEI) = DX*(I-1) + XORIGN
YI(NNODEI) = DY*(J-1) + YORIGN
END DO
END DO
WRITE(*,*) 'NUMBER OF DOMAIN ELEMENTS =',NEI
WRITE(*,*) 'NUMBER OF DOMAIN NODAL POINTS =',NNODEI
RETURN
END
C
C
SUBROUTINE SIZES ( ND,MXE,MXN,NE,NNODE,NODEX,X,Y,DS,
* BODYX,BODYY,CENTERX, CENTERY )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION NODEX(MXE,ND), X(MXN),Y(MXN)
C
C------ BODYX = SIZE OF OBJECT IN X COORDINATE
C------ BODYY = SIZE OF OBJECT IN Y COORDINATE
C
C------ (CENTERX,CENTERY) = CENTER OF OBJECT
C
C------ AVERAGE ELEMENT LENGTH
XMIN = X(1)
XMAX = X(1)
YMIN = Y(1)
YMAX = Y(1)
DO I = 2 , NNODE
XMIN = DMIN1 ( XMIN, X(I) )
XMAX = DMAX1 ( XMAX, X(I) )
YMIN = DMIN1 ( YMIN, Y(I) )
YMAX = DMAX1 ( YMAX, Y(I) )
END DO
BODYX = XMAX - XMIN
BODYY = YMAX - YMIN
CENTERX = ( XMAX + XMIN ) / 2.D0
CENTERY = ( YMAX + YMIN ) / 2.D0
DS = 0.D0
DO I = 1 , NE
DX = X(NODEX(I,2)) - X(NODEX(I,1))
DY = Y(NODEX(I,2)) - Y(NODEX(I,1))
DS = DS + DSQRT ( DX*DX + DY*DY )
END DO
DS = DS / DFLOAT(NE)
RETURN
END
C
C
SUBROUTINE WAVEHT ( OMEGA, MXNI,NNODEI,HI, XI, WAVEH )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION WAVEH(MXNI),XI(MXNI)
COMPLEX*16 HI(MXNI)
DO I = 1 , NNODEI
REALPT = DREAL ( HI(I) ) + DCOS( OMEGA*XI(I) )
FIMGPT = DIMAG ( HI(I) ) + DSIN( OMEGA*XI(I) )
WAVEH(I) = DSQRT ( REALPT**2 + FIMGPT**2 )
END DO
RETURN
END
C
C
SUBROUTINE INPUT (OMEGA,ND,MXE,MXN,NE,NNODE,NODEX,BV,X,Y,H,WAVEHB)
IMPLICIT REAL*8 ( A-H , O-Z )
COMPLEX*16 BV(MXE,ND),H(MXN)
DIMENSION X(MXE),Y(MXE),NODEX(MXE,ND),WAVEHB(MXN)
C=========>TAKING IN DATA FROM BEM1.DAT
OPEN ( 1, FILE='BEM1.DAT', STATUS='OLD' )
READ (1,*) OMEGA
READ (1,*) NE
DO IEL = 1 , NE
READ (1,*) I,(NODEX(I,J),J=1,ND)
END DO
READ(1,*) NNODE
DO I = 1 , NNODE
READ (1,*) NODE, X(NODE), Y(NODE)
END DO
CLOSE (1)
C=========> WAVEHEIGHT ON BOUNDARY NODES
OPEN ( 1, FILE='WHCOEFFI.CNT', STATUS='OLD' )
DO I = 1 , NNODE
READ (1,*) DUMMY, WAVEHB(I)
END DO
CLOSE (1)
C=========> TAKING IN BOUNDARY VALUES COMPUTED FROM ((HLM8LINQ.FOR))
OPEN (1,FILE='BOUNDARY.BIN',STATUS='UNKNOWN',FORM='UNFORMATTED')
DO I = 1 , NE
READ(1) BV(I,1), BV(I,2)
END DO
DO I = 1 , NNODE
READ(1) H(I)
END DO
CLOSE (1)
RETURN
END
C
C
SUBROUTINE INTE (PI,OMEGA,INTEPT,ND,XP,YP,C2,XE,YE,SAI,W,GE,FE)
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION XE(ND),YE(ND), SAI(INTEPT),W(INTEPT)
COMPLEX*16 GE(ND),FE(ND),C2, H0, H1
DO I = 1 , ND
GE(I) = DCMPLX(0.D0,0.D0)
FE(I) = DCMPLX(0.D0,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 )
C------- COMPUTATION OF SHAPE FUNCTIONS
SF1 = 0.5D0 * ( 1.D0 - SAI(IGAUSS) )
SF2 = 0.5D0 * ( 1.D0 + SAI(IGAUSS) )
C-------- COMPUTATION OF BESSEL FUNCTIONS --------
Z = R * OMEGA
CALL BESSEL ( PI, Z , H0, H1 )
C-------- INTEGRATION OF G(R)
GE(1) = GE(1) + H0 * W(IGAUSS) * SF1
GE(2) = GE(2) + H0 * W(IGAUSS) * SF2
C-------- INTEGRATION OF F(R)
RDOTN = (RX*DY-RY*DX) / DS
DRDN = RDOTN / R
FE(1) = FE(1) + OMEGA*H1 * DRDN * W(IGAUSS) * SF1
FE(2) = FE(2) + OMEGA*H1 * DRDN * W(IGAUSS) * SF2
END DO
GE(1) = C2 * DETJ * GE(1)
GE(2) = C2 * DETJ * GE(2)
FE(1) = -C2 * DETJ * FE(1)
FE(2) = -C2 * DETJ * FE(2)
RETURN
END
C
C
SUBROUTINE DOMAIN ( PI,OMEGA,INTEPT,ND,MXE,MXN,MXNI,C2,NNODEI,
* NE, GE,FE, SAI,W,XI,YI,NODEX,X,Y,H,BV,HI )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION NODEX(MXE,ND), SAI(INTEPT), W(INTEPT),
* X(MXN),Y(MXN),XE(ND),YE(ND)
DIMENSION XI(MXNI),YI(MXNI)
COMPLEX*16 GE(ND),FE(ND), BV(MXE,ND),H(MXN), SUM
COMPLEX*16 HI(MXNI), C2
C
WRITE(*,*) 'THIS PROCESS REQUIRES A LOT OF CPU TIME.'
C1 = - 1.D0/ ( 8.D0* DATAN( 1.D0) )
DO INSIDE = 1 , NNODEI
IF(INSIDE-INSIDE/100*100 .EQ. 0 ) THEN
IPERCENT= 100.- FLOAT(INSIDE)/FLOAT(NNODEI)*100.
WRITE(*,*)'REMAINING CALCULATION TIME = ',IPERCENT, '%'
END IF
XP = XI(INSIDE)
YP = YI(INSIDE)
SUM = DCMPLX(0.D0,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 INTELPS ( INTEPT,ND,XP,YP,C1,XE,YE,SAI,W, FE1,FE2 )
C = C + FE1 + FE2
END DO
IF ( C .GT. 0.9D0 ) THEN
DO IEL = 1 , NE
DO I = 1 , ND
XE(I) = X( NODEX(IEL,I) )
YE(I) = Y( NODEX(IEL,I) )
END DO
CALL INTE ( PI,OMEGA,INTEPT,ND,XP,YP,C2,XE,YE,SAI,W,GE,FE )
DO J = 1 , ND
SUM = SUM + FE(J)*H(NODEX(IEL,J)) - GE(J)*BV(IEL,J)
END DO
END DO
END IF
HI(INSIDE) = SUM/C
END DO
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'
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
IF ( N .EQ. 5 ) THEN
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
SAI(4) = -SAI(2)
SAI(5) = -SAI(1)
W(4) = W(2)
W(5) = W(1)
RETURN
END IF
IF ( N .EQ. 6 ) THEN
SAI(1) = 0.23861918608320D0
SAI(2) = 0.66120938646626D0
SAI(3) = 0.93246951420315D0
W(1) = 0.46791393457269D0
W(2) = 0.36076157304814D0
W(3) = 0.17132449237917D0
SAI(4) = -SAI(1)
SAI(5) = -SAI(2)
SAI(6) = -SAI(3)
W(4) = W(1)
W(5) = W(2)
W(6) = W(3)
END IF
RETURN
END
C
C
SUBROUTINE BESSEL ( PI, X , H0, H1 )
IMPLICIT REAL*8 ( A-H , O-Z )
COMPLEX*16 H0, H1
PARAMETER ( EULER = 0.5772156649015D0, ALLOWANC=1.D-16 )
C--------- BJ0, BJ1
X22 = (X/2.D0) **2
PRODUCT = 1.D0
SIGNCH = 1.D0
FACT = 1.D0
SUM0 = 0.D0
SUM1 = 0.D0
M = 1
CHECK = 1.D0
DO WHILE ( DABS(CHECK) .GT. ALLOWANC )
FACT = FACT * DFLOAT(M)
PRODUCT = PRODUCT * X22
SIGNCH = -SIGNCH
DELTA0 = SIGNCH*PRODUCT/(FACT*FACT)
SUM0 = SUM0 + DELTA0
IF ( SUM0 .EQ. 0.D0 ) THEN
CHECK = DELTA0
ELSE
CHECK = DELTA0 / SUM0
END IF
SUM1 = SUM1 + SIGNCH*PRODUCT/(FACT*FACT*DFLOAT(M+1))
M = M + 1
IF ( M .GT. 1000 ) STOP
END DO
BJ0 = 1.D0 + SUM0
BJ1 = (X/2.D0) * ( 1.D0 + SUM1 )
C
C-------- BN0, BN1
PRODUCT = 1.D0
SIGNCH = 1.D0
FACT = 1.D0
SUM0 = 0.D0
SUM1 = 0.D0
M = 1
CHECK = 1.D0
DO WHILE ( DABS(CHECK) .GT. ALLOWANC )
FACT = FACT * DFLOAT(M)
PRODUCT = PRODUCT * X22
SIGNCH = -SIGNCH
C..............SUM OF 1/1+1/2+1/3+.......+1/M
SUMM = 0.D0
DO I = 1 , M
SUMM = SUMM + 1.D0 / DFLOAT(I)
END DO
DELTA0 = SIGNCH*PRODUCT/(FACT*FACT) * 2.D0*SUMM
SUM0 = SUM0 + DELTA0
IF ( SUM0 .EQ. 0.D0 ) THEN
CHECK = DELTA0
ELSE
CHECK = DELTA0 / SUM0
END IF
SUM1 = SUM1 + SIGNCH*PRODUCT/(FACT*FACT*DFLOAT(M+1)) *
* (2.D0*SUMM+1.D0/DFLOAT(M+1))
M = M + 1
IF ( M .GT. 1000 ) STOP
END DO
BN0 = (2.D0/PI)*BJ0*(EULER+DLOG(X/2.D0)) - 1.D0/PI*SUM0
BN1 = (2.D0/PI)*BJ1*(EULER+DLOG(X/2.D0)) - X/(2.D0*PI) -
* 1.D0/PI*(X/2.D0)*SUM1 - 2.D0/(PI*X)
C----------- ASSEMBLE THEM INTO COMPLEX NUMBERS
H0 = DCMPLX(BJ0,BN0)
H1 = DCMPLX(BJ1,BN1)
RETURN
END
C
C
SUBROUTINE INTELPS (INTEPT,ND,XP,YP,C1, XE, YE, SAI,W,FE1,FE2 )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION XE(ND),YE(ND), SAI(INTEPT),W(INTEPT)
FE1 = 0.D0
FE2 = 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--------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 F(R)
TEMP = (RX*DY-RY*DX) / (R*R) * W(IGAUSS)
FE1 = FE1 + TEMP * SF1
FE2 = FE2 + TEMP * SF2
END DO
FE1 = -C1 * DETJ /DS * FE1
FE2 = -C1 * DETJ /DS * FE2
RETURN
END