PROGRAM GRA4FILE
C=======================================================================
C   FINITE ELEMENT GRAPHICS FOR TWO-DIMENSIONAL LAPLACE SOLUTION
C            ELEMENT TYPE: FOUR-NODED ISOPARAMETRIC ELEMENT
C              INPUT FILE: FEM04INPUT.DAT
C             OUTPUT FILE: TEMPERAT.OUT
C         ORIGINAL CODE: NSEQ8GRA.FOR, EIJI FUKUMORI, JULY 1985
C 25JAN. 2013
C=======================================================================
      INCLUDE 'PARAM.DAT'
      IMPLICIT REAL*8 ( A-H , O-Z )
C                                ARRAYS
      CHARACTER*14 INPFILE
      DIMENSION SK(ND,ND),BX(2,ND), S(ND)
      DIMENSION NODEX(MXE,ND), XCOORD(2,MXN),U(MXN)
C=======================================================================
      WRITE (*,*)' READING IN DATA FILES'
      INPFILE = 'FEM04INPUT.DAT'
      CALL INPUT ( INPFILE,ND,MXE,MXN,NE,NNODE,EXX,EXY,EYY,NODEX,
     * XCOORD,U )
      WRITE (*,*)'PROJECT NAME =======>', INPFILE
C=======================================================================
      NSTEP = 20
      CALL TEMPA (MXE,MXN,ND,NE,NNODE,NODEX,XCOORD,S,BX,U,NSTEP )
C=======================================================================
      STOP 'NORMAL TERMINATION'
      END
C
C
      SUBROUTINE INPUT ( INPFILE,ND,MXE,MXN,NE,NNODE,EXX,EXY,EYY,
     *   NODEX, XCOORD, U )
      IMPLICIT REAL*8 ( A-H , O-Z )
      DIMENSION NODEX(MXE,ND),XCOORD(2,MXN), U(MXN)
      CHARACTER INPFILE*14
      OPEN ( 1, FILE = INPFILE, STATUS = 'UNKNOWN' )
      READ (1,*) EXX, EXY, EYY
      READ (1,*) NE
      IF ( NE .GT. MXE ) STOP 'NE > MXE'
      READ (1,*) (IEL,(NODEX(IEL,J),J=1,ND), I=1,NE)
      READ (1,*) NNODE
      IF ( NNODE .GT. MXN ) STOP 'NNODE > MXN'
      READ (1,*) (NODE,XCOORD(1,NODE),XCOORD(2,NODE),J=1,NNODE)
C=======================================================================
C========> FILENAME XXXXXXX.BIN
      OPEN (1,FILE="SOLUTION.BIN",STATUS='UNKNOWN',FORM='UNFORMATTED')
      READ (1) ( U(I) , I = 1 , NNODE )
      CLOSE (1)
      RETURN
      END
C
C
      SUBROUTINE PLTSAI ( ND, DS, NSTEP, START, CRD, S )
      IMPLICIT REAL*8 ( A-H , O-Z )
      DIMENSION CRD(2,ND), S(ND)
      IF ( NSTEP .EQ. 0 ) RETURN
      DO LEVEL = 1 , NSTEP
      SXY = START + (LEVEL-1) * DS
      K = 1
      DO I = 1 , ND
      J = I + 1
      IF ( I .EQ. ND ) J = 1
      IF ( (S(I)-SXY)*(S(J)-SXY) .LT. 0 ) THEN
           T=(SXY-S(I))/(S(J)-S(I))
           X0=CRD(1,J)*T+(1.D0-T)*CRD(1,I)
           Y0=CRD(2,J)*T+(1.D0-T)*CRD(2,I)
           IF ( K .GT. 0 ) THEN
                CALL XMOVE ( X0, Y0 )
           ELSE
                CALL XDRAW ( X0, Y0 )
           END IF
           K = -K
      END IF
      END DO
      END DO
      RETURN
      END
C
C
      SUBROUTINE CONTOUR (MXE,MXN,ND,NE,NNODE,NODEX,XCOORD,S,B,P,NSTEP)
      IMPLICIT REAL*8 ( A-H , O-Z )
      DIMENSION S(ND),NODEX(MXE,ND),XCOORD(2,MXN),B(2,ND),P(MXN)
      CALL MINMAX ( MXN, NNODE, P, PPMIN, PPMAX )
      IF ( PPMAX .EQ. PPMIN ) RETURN
      NSTEP = NSTEP/2*2 + 1
      DS = ( PPMAX - PPMIN ) / NSTEP
      I = PPMIN/DS
      PPMIN= I*DS
C------------
      CALL BOUND (MXE,MXN,ND,NE,NODEX,XCOORD)
      DO IEL = 1 , NE
      DO I   = 1 , ND
      B(1,I) = XCOORD(1,NODEX(IEL,I))
      B(2,I) = XCOORD(2,NODEX(IEL,I))
      S(I)   =        P(NODEX(IEL,I))
      END DO
      CALL PLTSAI ( ND, DS, NSTEP, PPMIN, B, S )
      END DO
      RETURN
      END
C
C
      SUBROUTINE TEMPA (MXE,MXN,ND,NE,NNODE,NODEX,XCOORD,S,B,P,NSTEP )
      IMPLICIT REAL*8 ( A-H , O-Z )
      DIMENSION NODEX(MXE,ND),XCOORD(2,MXN),P(MXN),S(ND),B(2,ND)
      CHARACTER FILENAME*12
      FILENAME = "TEMPERAT.OUT"
      CALL PLTLGO ( FILENAME )
      CALL CONTOUR (MXE,MXN,ND,NE,NNODE,NODEX,XCOORD,S,B,P,NSTEP)
      CALL PLTEXT
      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 = DMIN1 ( QMIN , Q(I) )
      QMAX = DMAX1 ( QMAX , Q(I) )
      END DO
      RETURN
      END
C
C
      SUBROUTINE BOUND (MXE,MXN,ND,NE,NODEX,XCOORD)
      IMPLICIT REAL*8 ( A-H , O-Z )
      DIMENSION XCOORD(2,MXN),NODEX(MXE,ND)
      LOGICAL LINE
C
      DO IEL = 1 , NE
      DO I = 1 , ND
      LINE = .TRUE.
      J = I + 1
      IF ( I .EQ. ND ) J = 1
C
      DO JEL = 1 , NE
      IF ( IEL .NE. JEL ) THEN
      DO II = 1 , ND
      JJ = II + 1
      IF ( II .EQ. ND ) JJ = 1
      IF ( NODEX(IEL,I).EQ.NODEX(JEL,JJ) .AND.
     *     NODEX(IEL,J).EQ.NODEX(JEL,II) ) THEN
      LINE = .FALSE.
      EXIT
      END IF
      END DO
      END IF
      IF ( .NOT. LINE ) EXIT
      END DO
      IF ( LINE ) THEN
      CALL XMOVE ( XCOORD(1,NODEX(IEL,I)) , XCOORD(2,NODEX(IEL,I)) )
      CALL XDRAW ( XCOORD(1,NODEX(IEL,J)) , XCOORD(2,NODEX(IEL,J)) )
      END IF
C
      END DO
      END DO
      RETURN
      END
C
C
C======================== GRAPHICS RUTINES =============================
      SUBROUTINE PLTLGO ( FILENAME )
      IMPLICIT REAL*8 ( A-H , O-Z )
      CHARACTER FILENAME*12
      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