PROGRAM LAGFILE
C=======================================================================
C FINITE ELEMENT GRAPHICS FOR TWO-DIMENSIONAL POTENTIAL PROBLEMS
C ELEMENT TYPE: FOUR-NODED ISOPARAMETRIC ELEMENT
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)
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,IBND,ITYPE, BV )
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 )
IF ( IITYPE .EQ. 2 )
* CALL POTENTL (MXE,MXN,ND,NE,NNODE,NODEX,XCOORD,YCOORD,POTEN,
* S,BX,IB)
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,IBND,ITYPE, BV )
DIMENSION NODEX(MXE,ND),XCOORD(MXN),YCOORD(MXN),IBND(MXB),
* POTEN(MXN),ITYPE(MXB),BV(MXB)
LOGICAL YES
C========> FILE INQUIRE
INQUIRE ( FILE='HEL4DATA.QQQ', EXIST=YES )
IF ( .NOT.YES ) STOP 'FILE DOES NOT EXIST'
OPEN ( 1, FILE='HEL4DATA.QQQ', STATUS='OLD' )
READ (1,*) ALPHA
C========> ELEMENT
READ (1,*) NE
DO I = 1 , NE
READ (1,*) IEL, ( NODEX(IEL,J), J = 1 , ND )
END DO
C========> COORDINATES
READ (1,*) NNODE
DO I = 1 , NNODE
READ (1,*) NODE, XCOORD(NODE) , YCOORD(NODE)
END DO
C========> BOUNDARY CONDITIONS
READ(1,*) NB
DO I = 1 , NB
READ (1,*) IBND(I) , ITYPE(I), BV(I)
END DO
CLOSE (1)
C========> POTENTIAL VALUE FILENAME=SOLUTION.BIN
OPEN (1,FILE='SOLUTION.BIN',STATUS='UNKNOWN',FORM='UNFORMATTED')
READ (1) ( POTEN(I) , I = 1 , NNODE )
CLOSE (1)
RETURN
END
C
C
SUBROUTINE POTENTL ( MXE,MXN,ND,NE,NNODE,NODEX, XCOORD,YCOORD,
* POTEN,S,B,IB )
DIMENSION S(ND),NODEX(MXE,ND),XCOORD(MXN),YCOORD(MXN),IB(MXE,ND),
* B(2,ND),POTEN(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 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 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)
DIMENSION NODEX(MXE,ND),XCOORD(MXN),YCOORD(MXN),IB(MXE,ND)
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
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