PROGRAM WAVEPLOT
C=======================================================================
C FINITE ELEMENT GRAPHICS FOR
C TWO-DIMENSIONAL POTENTIAL PROBLEMS
C ELEMENT TYPE: FOUR-NODED ELEMENT FOR DOMAIN
C TWO-NODED ELEMENT FOR BOUNDARY
C CREATE FILES CONTAINING XY-COORDINTE VALUES
C ORIGINAL CODE: NSEQ8GRA.FOR, EIJI FUKUMORI, JULY 1985
C=======================================================================
IMPLICIT REAL*8 ( A-H , O-Z )
PARAMETER ( MXE=30000,MXN=33000,MXB=10000,ND=4)
C ARRAYS
DIMENSION NODEX(MXE,ND),ISEG(MXE,ND),IB(MXE,ND)
DIMENSION IBND(MXB), BV(MXB), ITYPE(MXB)
DIMENSION XCOORD(MXN),YCOORD(MXN),POTEN(MXN)
DIMENSION S(ND), BX(2,ND), NODEXB(MXE,ND)
DIMENSION XB(MXN), YB(MXN), PP(MXN)
C=======================================================================
WRITE(*,*)' POTENTIAL GRAPHICS PROGRAM'
WRITE (*,*)' READING IN DATA FILES'
C
CALL INPUT ( MXE,MXN,MXB,ND,NE,NNODE,NB,NODEX,XCOORD,YCOORD,
* POTEN,XB,YB,PP,NEB,NNODEB,NODEXB )
C=======================================================================
CALL PLTEL4 (ND,MXE,MXN,NODEX,NE,XCOORD,YCOORD,NNODE,IB,
* NEB,NODEXB,XB,YB )
C=====================================================================
WRITE(*,*)' TYPE IN NUMBER OF CONTOUR LINES(N=20 ABOUT RIGHT)'
NSTEP = 20
WRITE(*,*)' N=', NSTEP
C=====================================================================
CALL POTENTL (MXE,MXN,ND,NE,NNODE,NODEX,XCOORD,YCOORD,POTEN,
* S,BX,IB, NEB,NODEXB,XB,YB, NSTEP )
STOP 'TERMINATION'
END
C
C
SUBROUTINE INPUT ( MXE,MXN,MXB,ND,NE,NNODE,NB,NODEX,XCOORD,
* YCOORD,POTEN,XB,YB,PP,NEB,NNODEB,NODEXB )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION NODEX(MXE,ND),XCOORD(MXN),YCOORD(MXN),IBND(MXB),
* POTEN(MXN),ITYPE(MXB),BV(MXB),NODEXB(MXE,ND)
DIMENSION XB(MXN),YB(MXN),PP(MXN)
C
OPEN (1,FILE='DOMAIN.BIN',STATUS='UNKNOWN',FORM='UNFORMATTED')
READ(1) NE
READ(1) NNODE
DO I = 1 , NE
READ(1) (NODEX(I,J),J=1,ND)
END DO
READ(1) ( XB(I) , I=1 , NNODE )
READ(1) ( YB(I) , I=1 , NNODE )
READ(1) ( PP(I) , I=1 , NNODE )
CLOSE (1)
DO I = 1 , NNODE
XCOORD(I) = XB(I)
YCOORD(I) = YB(I)
POTEN (I) = PP(I)
END DO
WRITE(*,*) 'NUMBER OF DOMAIN ELEMENTS =',NE
WRITE(*,*) 'NUMBER OF DOMAIN NODAL POINTS =',NNODE
C
OPEN ( 1, FILE='BEM1.DAT', STATUS='UNKNOWN' )
READ (1,*) OMEGA
READ (1,*) NEB
DO IEL = 1 , NEB
READ (1,*) I,(NODEXB(I,J),J=1,2)
END DO
READ(1,*) NNODEB
DO I = 1 , NNODEB
READ (1,*) NODE, XB(NODE), YB(NODE)
END DO
CLOSE (1)
C
RETURN
END
C
C
SUBROUTINE POTENTL ( MXE,MXN,ND,NE,NNODE,NODEX, XCOORD,YCOORD,
* POTEN,S,B,IB,NEB,NODEXB,XB,YB, NSTEP )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION S(ND),NODEX(MXE,ND),XCOORD(MXN),YCOORD(MXN),IB(MXE,ND),
* B(2,ND),POTEN(MXN),NODEXB(MXE,ND), XB(MXN),YB(MXN)
CHARACTER FILENAME*12
COMMON / PLOTTING / FILENAME
COMMON / PORT / TMIN, TMAX / INCREM / DT
FILENAME = 'POTENTAL.OUT'
CALL MINMAX ( MXN, NNODE, POTEN, TMIN, TMAX )
CALL CONTOUR (MXE,MXN,ND,NE,NNODE,NODEX,XCOORD,YCOORD,
* S,B,IB,POTEN, NSTEP)
CALL BOUNDB ( MXE,MXN,ND,NEB,NODEXB,XB,YB )
CALL PLTEXT
RETURN
END
C
C
SUBROUTINE BOUNDB ( MXE,MXN,ND,NEB,NODEXB,XB,YB )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION NODEXB(MXE,ND), XB(MXN),YB(MXN)
DO I = 1 , NEB
X1 = XB(NODEXB(I,1))
Y1 = YB(NODEXB(I,1))
X2 = XB(NODEXB(I,2))
Y2 = YB(NODEXB(I,2))
CALL XMOVE ( X1 , Y1 )
CALL XDRAW ( X2 , Y2 )
END DO
RETURN
END
C
C
SUBROUTINE CONTOUR (MXE,MXN,ND,NE,NNODE,NODEX,XCOORD,YCOORD,
* S, B, IB, P, NSTEP )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION NODEX(MXE,ND), XCOORD(MXN),YCOORD(MXN), P(MXN),
* S(ND), B(2,ND), IB(MXE,ND)
COMMON / PORT / PPMIN, PPMAX, / INCREM / DS
IF ( PPMAX .EQ. PPMIN ) RETURN
NSTEP = NSTEP/2*2
DS = ( PPMAX - PPMIN ) / NSTEP
C
IF ( (PPMIN .LT. 0.) .AND. (PPMAX.GT.0.) ) THEN
DS = AMAX1(PPMAX, -PPMIN)/NSTEP
END IF
WRITE(*,*) 'NSTEP =',NSTEP, ' DS = ',DS
WRITE(*,*)' PPMAX=',PPMAX,' PPMIN=',PPMIN
C
CALL PLTLGO
IC = 0
CALL BOUND ( MXE,MXN,ND,NE,NNODE,NODEX,XCOORD,YCOORD,IB,IC )
DO IEL = 1 , NE
DO I = 1 , ND
B(1,I) = XCOORD(NODEX(IEL,I))
B(2,I) = YCOORD(NODEX(IEL,I))
S(I) = P(NODEX(IEL,I))
END DO
IF ( (PPMIN .LT. 0.) .AND. (PPMAX.GT.0.) ) THEN
CALL PLTSAI ( DS, NSTEP/2, -NSTEP*DS/2, B, S )
C CALL PLTSAI ( DS, 1, 0., B, S )
CALL PLTSAI ( DS, NSTEP/2, DS, B, S )
ELSE
CALL PLTSAI ( DS, NSTEP, PPMIN, B, S )
END IF
END DO
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 = AMAX1 ( SS(1), SS(2), SS(3), SS(4) )
SMIN = AMIN1 ( SS(1), SS(2), SS(3), SS(4) )
IF ( START .GT. SMAX ) RETURN
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 BOUND (MXE,MXN,ND,NE,NNODE,NODEX,XCOORD,YCOORD,IB,IC)
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION XCOORD(MXN),YCOORD(MXN),NODEX(MXE,ND),IB(MXE,ND)
C---------- IF IC=0, BOUNDARY LINE ONLY
C---------- IF IC=1, ELEMENTS
C--------- INITIALIZATION
IF ( IC .EQ. 1 ) THEN
DO IEL = 1 , NE
DO J = 1 , ND
IS = J
IE = J + 1
IF ( IS .EQ. ND ) IE = 1
NDIS = NODEX(IEL,IS)
NDIE = NODEX(IEL,IE)
CALL XMOVE ( XCOORD(NDIS) , YCOORD(NDIS) )
CALL XDRAW ( XCOORD(NDIE) , YCOORD(NDIE) )
END DO
END DO
RETURN
END IF
C=============== FOR A CASE OF IC=0 ==================
DO I = 1 , NE
DO J = 1 , ND
IB(I,J) = 0
END DO
END DO
C
DO IEL = 1 , NE-1
DO ISEG = 1 , ND
IS = ISEG
IE = ISEG+1
IF ( IS .EQ. ND ) IE = 1
NDI1 = NODEX(IEL,IS)
NDI2 = NODEX(IEL,IE)
IF ( IB(IEL,ISEG) .EQ. 0 ) THEN
DO JEL = IEL+1 , NE
DO JSEG = 1 , ND
JS = JSEG
JE = JSEG+1
IF ( JS .EQ. ND ) JE = 1
NDJ1 = NODEX(JEL,JS)
NDJ2 = NODEX(JEL,JE)
IF ( (NDI1 .EQ. NDJ2) .AND. (NDI2 .EQ. NDJ1) ) THEN
IB(JEL,JSEG) = 1
IB(IEL,ISEG) = 1
END IF
END DO
END DO
END IF
END DO
END DO
C
DO IEL = 1 , NE
DO ISEG = 1 , ND
IF ( IB(IEL,ISEG) .EQ. 0 ) THEN
IS = ISEG
IE = ISEG+1
IF ( IS .EQ. ND ) IE = 1
IS = NODEX(IEL,IS)
IE = NODEX(IEL,IE)
CALL XMOVE ( XCOORD(IS) , YCOORD(IS) )
CALL XDRAW ( XCOORD(IE) , YCOORD(IE) )
END IF
END DO
END DO
RETURN
END
C
C
SUBROUTINE MINMAX ( MXN, NNODE, Q, QMIN, QMAX )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION Q(MXN)
QMIN = Q(1)
QMAX = Q(1)
DO I = 1 , NNODE
QMIN = AMIN1 ( QMIN , Q(I) )
QMAX = AMAX1 ( QMAX , Q(I) )
END DO
RETURN
END
C
C
SUBROUTINE PLTEL4 (ND,MXE,MXN,NODEX,NE,XCOORD,YCOORD,NNODE,IB,
* NEB,NODEXB,XB,YB )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION NODEX(MXE,ND),XCOORD(MXN),YCOORD(MXN),IB(MXE,ND),
* NODEXB(MXE,ND), XB(MXN),YB(MXN)
CHARACTER FILENAME*12
COMMON / PLOTTING / FILENAME
FILENAME = "ELEMENT.OUT"
CALL PLTLGO
CALL BOUNDB ( MXE,MXN,ND,NEB,NODEXB,XB,YB )
IC = 1
CALL BOUND ( MXE,MXN,ND,NE,NNODE,NODEX,XCOORD,YCOORD,IB,IC )
CALL PLTEXT
RETURN
END
C
C
SUBROUTINE SEGCHK ( ND,MXE,NODEX,NE,ISEG,II,JJ )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION NODEX(MXE,ND), ISEG(MXE,ND)
DO IEL = 1 , NE
DO I = 1 , ND
IF ( NODEX(IEL,I) .EQ.JJ ) THEN
J = I + 1
IF ( I .EQ. ND ) J = 1
IF ( NODEX(IEL,J).EQ.II ) THEN
ISEG (IEL,I) = 0
RETURN
ENDIF
ENDIF
END DO
END DO
RETURN
END
C
C======================== GRAPHICS RUTINES =============================
SUBROUTINE PLTLGO
IMPLICIT REAL*8 ( A-H , O-Z )
CHARACTER FILENAME*12
COMMON / PLOTTING / FILENAME
OPEN ( 1, FILE=FILENAME, STATUS='UNKNOWN' )
RETURN
END
C
C
SUBROUTINE PLTEXT
IMPLICIT REAL*8 ( A-H , O-Z )
CLOSE (1)
RETURN
END
C
C
SUBROUTINE XMOVE ( X , Y )
IMPLICIT REAL*8 ( A-H , O-Z )
WRITE(1,*) X, Y
RETURN
END
C
C
SUBROUTINE XDRAW ( X , Y )
IMPLICIT REAL*8 ( A-H , O-Z )
WRITE(1,*) X, Y
WRITE(1,*)
RETURN
END