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=======================================================================
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)
REAL*8 XB(MXN), YB(MXN), PP(MXN)
COMMON / DOMAIN / XMIN, XMAX, YMIN , YMAX
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 )
CALL MINMAX ( MXN, NNODE, XCOORD, XMIN , XMAX )
CALL MINMAX ( MXN, NNODE, YCOORD, YMIN , YMAX )
WRITE (*,*)' MAX & MIN IN X-COORDINATE: ', XMAX, XMIN
WRITE (*,*)' MAX & MIN IN Y-COORDINATE: ', YMAX, YMIN
IF ( XMAX-XMIN .LE. 0. ) STOP' OBJECT HAS NO LENGTH IN XCOORD.'
IF ( YMAX-YMIN .LE. 0. ) STOP' OBJECT HAS NO LENGTH IN YCOORD.'
C=======================================================================
NMENU = 2
10 CALL MENU ( IITYPE )
IF ( IITYPE .EQ. 1 )
* CALL PLTEL4 (ND,MXE,MXN,NODEX,NE,XCOORD,YCOORD,NNODE,IB,
* NEB,NODEXB,XB,YB )
IF ( IITYPE .EQ. 2 )
* CALL POTENTL (MXE,MXN,ND,NE,NNODE,NODEX,XCOORD,YCOORD,POTEN,
* S,BX,IB, NEB,NODEXB,XB,YB )
IF (IITYPE.GT.NMENU ) STOP 'TERMINATION'
IF (IITYPE.LE.0) STOP 'TERMINATION'
GO TO 10
END
C
C
SUBROUTINE MENU ( ITYPE )
WRITE (*,*)'----------------------------------------------------'
WRITE (*,*)' ID # MENU'
WRITE (*,*)'===================================================='
WRITE (*,*)' 0 OR LESS TERMINATION'
WRITE (*,*)' 1 FINITE ELEMENT DISCRETIZATION'
WRITE (*,*)' 2 POTENTIAL CONTOUR'
WRITE (*,*)' 3 OR MORE TERMINATION'
WRITE (*,*)'----------------------------------------------------'
WRITE (*,*)' TYPE IN ID # = '
READ (*,*) ITYPE
RETURN
END
C
C
SUBROUTINE INPUT ( MXE,MXN,MXB,ND,NE,NNODE,NB,NODEX,XCOORD,
* YCOORD,POTEN,XB,YB,PP,NEB,NNODEB,NODEXB )
DIMENSION NODEX(MXE,ND),XCOORD(MXN),YCOORD(MXN),IBND(MXB),
* POTEN(MXN),ITYPE(MXB),BV(MXB),NODEXB(MXE,ND)
REAL*8 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 )
DIMENSION S(ND),NODEX(MXE,ND),XCOORD(MXN),YCOORD(MXN),IB(MXE,ND),
* B(2,ND),POTEN(MXN),NODEXB(MXE,ND)
REAL*8 XB(MXN),YB(MXN)
CHARACTER FILENAME*8
COMMON / PLOTTING / FILENAME
COMMON / PORT / TMIN, TMAX / INCREM / DT
FILENAME = 'POTENTAL'
CALL MINMAX ( MXN, NNODE, POTEN, TMIN, TMAX )
CALL CONTOUR (MXE,MXN,ND,NE,NNODE,NODEX,XCOORD,YCOORD,
* S,B,IB,POTEN)
TREF = 0.
TIME = 0.
WRITE(*,200) TREF,TIME,TMIN, TMAX, DT
CALL BOUNDB ( MXE,MXN,ND,NEB,NODEXB,XB,YB )
CALL PLTEXT
200 FORMAT ( 1X,'POTENTIAL FIELD:'/ 1X,'BLUE=TREF'/
* 1X, 'RED=BELOW TREF'/ 1X, 'GREEN=ABOVE TREF'/1X,'TREF=',G10.3,
* 19(/),1X, 'CURRENT VALUES:'/1X,'TIME=',G10.3 / 1X, 'MIN=',G10.3/
* 1X, 'MAX=',G10.3 / 1X,'DS=',G10.3 )
RETURN
END
C
C
SUBROUTINE BOUNDB ( MXE,MXN,ND,NEB,NODEXB,XB,YB )
DIMENSION NODEXB(MXE,ND)
REAL*8 XB(MXN),YB(MXN)
REAL*4 X1, X2, Y1, Y2
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 )
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
WRITE(*,*)' TYPE IN NUMBER OF CONTOUR LINES(N=20 ABOUT RIGHT)'
WRITE(*,*)' N=0 FOR QUIT, MAX N = 99'
WRITE(*,*)' N='
READ (*,*) NSTEP
IF ( NSTEP .LE. 0 ) RETURN
IF ( NSTEP .GT. 99 ) 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 )
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)
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
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
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)
DO JEL = IEL+1 , NE
DO JSEG = 1 , ND
IF ( IB(JEL,JSEG) .EQ. 0 ) THEN
JS = JSEG
JE = JSEG+1
IF ( JS .EQ. ND ) JE = 1
JS = NODEX(JEL,JS)
JE = NODEX(JEL,JE)
IF ( (IS .EQ. JE) .AND. (IE .EQ. JS) ) THEN
IF ( IC .EQ. 0 ) IB(IEL,ISEG) = 1
IB(JEL,JSEG) = 1
END IF
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 )
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 )
DIMENSION NODEX(MXE,ND),XCOORD(MXN),YCOORD(MXN),IB(MXE,ND),
* NODEXB(MXE,ND)
REAL*8 XB(MXN),YB(MXN)
CHARACTER FILENAME*8
COMMON / PLOTTING / FILENAME
COMMON / DOMAIN / XMIN, XMAX, YMIN, YMAX
FILENAME = "ELEMENT"
CALL PLTLGO
WRITE (*,200) TIME, NE, NNODE, XMIN, XMAX, YMIN, YMAX
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
200 FORMAT ( 1X,'FINITE ELEMENT'/ 1X,'DISCRETIZATION'/
* 1X,'TIME=',G10.3 /1X,'NE=', I5 / 1X,'NNODE=', I5, 20(/),
* 1X,'XMIN=',G10.3 / 1X,'XMAX=',G10.3 /
* 1X,'YMIN=',G10.3 / 1X,'YMAX=',G10.3 )
RETURN
END
C
C
SUBROUTINE SEGCHK ( ND,MXE,NODEX,NE,ISEG,II,JJ )
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
CHARACTER FILENAME*8
COMMON / DOMAIN / XMIN , XMAX , YMIN , YMAX
COMMON / PENDOMAIN / IXMIN , IXMAX , IYMIN , IYMAX
COMMON / PLOTTING / FILENAME, RMAGNI, DIGITS
DIGITS = 99999.
OPEN ( 1, FILE=FILENAME, STATUS='UNKNOWN' )
RMAGNI = AMAX1 ( XMAX-XMIN, YMAX-YMIN )
IXMIN = DIGITS
IYMIN = DIGITS
IXMAX = -DIGITS
IYMAX = -DIGITS
RETURN
END
C
C
SUBROUTINE PLTEXT
COMMON / PENDOMAIN / IXMIN , IXMAX , IYMIN , IYMAX
WRITE(1,'(4I7)') IXMIN , IXMAX , IYMIN , IYMAX
CLOSE (1)
RETURN
END
C
C
SUBROUTINE XMOVE ( X , Y )
CHARACTER FILENAME*8
COMMON / DOMAIN / XMIN , XMAX , YMIN , YMAX
COMMON / PENDOMAIN / IXMIN , IXMAX , IYMIN , IYMAX
COMMON / PLOTTING / FILENAME, RMAGNI, DIGITS
COMMON / PENMOVE /IXMOVE, IYMOVE
IXMOVE = (X-XMIN)/RMAGNI*DIGITS
IYMOVE = (Y-YMIN)/RMAGNI*DIGITS
IXMIN = MIN0 ( IXMIN, IXMOVE )
IXMAX = MAX0 ( IXMAX, IXMOVE )
IYMIN = MIN0 ( IYMIN, IYMOVE )
IYMAX = MAX0 ( IYMAX, IYMOVE )
RETURN
END
C
C
SUBROUTINE XDRAW ( X , Y )
CHARACTER FILENAME*8
COMMON / PENDOMAIN / IXMIN , IXMAX , IYMIN , IYMAX
COMMON / DOMAIN / XMIN , XMAX , YMIN , YMAX
COMMON / PLOTTING / FILENAME, RMAGNI, DIGITS
COMMON / PENMOVE /IXMOVE, IYMOVE
IX = (X-XMIN)/RMAGNI*DIGITS
IY = (Y-YMIN)/RMAGNI*DIGITS
WRITE(1,'(4I7)') IXMOVE, IYMOVE, IX, IY
IXMIN = MIN0 ( IXMIN, IX )
IXMAX = MAX0 ( IXMAX, IX )
IYMIN = MIN0 ( IYMIN, IY )
IYMAX = MAX0 ( IYMAX, IY )
RETURN
END