PROGRAM NSG4FILE
C=======================================================================
C   FINITE ELEMENT GRAPHICS FOR TWO-DIMENSIONAL LAPLACE SOLUTIONS
C            ELEMENT TYPE: THREE-NODED LINEAR ELEMENT
C             CREATE FILE: TEMPERAT.OUT
C         ORIGINAL CODE: FEM3GRA.FOR, EIJI FUKUMORI, JULY 1985
C 13 FEB. 2013
C=======================================================================
      INCLUDE 'PARAM.DAT'
      IMPLICIT REAL*8 ( A-H , O-Z )
C                                ARRAYS
      DIMENSION SK(ND,ND),X(ND),Y(ND),BX(2,ND), S(ND)
      DIMENSION NODEX(MXE,ND)
      DIMENSION XCOORD(MXN), YCOORD(MXN),U(MXN)
C=======================================================================
      WRITE (*,*)' READING IN DATA FILES'
      CALL INPUT ( MXE,MXN,ND,NE,NNODE,NODEX,XCOORD,YCOORD,U )
C=======================================================================
      WRITE (*,*) 'CALL TEMPA'
      NSTEP = 25
      CALL TEMPA (MXE,MXN,ND,NE,NNODE,NODEX,XCOORD,YCOORD,S,BX,U,NSTEP)
C=======================================================================
      STOP 'NORMAL TERMINATION'
      END
C
C
      SUBROUTINE INPUT ( MXE,MXN,ND,NE,NNODE,NODEX,XCOORD,YCOORD,U )
      IMPLICIT REAL*8 ( A-H , O-Z )
      DIMENSION NODEX(MXE,ND),XCOORD(MXN),YCOORD(MXN),U(MXN)
      CHARACTER INPFILE*14
      INPFILE = 'FEM03INPUT.DAT'
      OPEN ( 1, FILE = INPFILE, STATUS = 'UNKNOWN' )
      READ (1,*) EXX, EXY, EYY
C========> ELEMENT
      READ (1,*) NE
                    IF ( NE    .GT. MXE) STOP'ERROR #1'
                    IF ( NE    .LE. 0  ) STOP'NE=0'
      DO I = 1 , NE
      READ (1,*) IEL, ( NODEX(IEL,J), J = 1 , ND )
      END DO
C========> NODE
      READ (1,*) NNODE
                    IF ( NNODE .GT. MXN) STOP'ERROR #2'
                    IF ( NNODE .LE. 0  ) STOP'NNODE=0'
      DO I = 1 , NNODE
      READ (1,*) NODE, XCOORD(NODE) , YCOORD(NODE)
      END DO
      CLOSE (1)
C========> FILENAME XXXXXXX.BIN
      OPEN (1,FILE="SOLUTION.BIN",STATUS='UNKNOWN',FORM='UNFORMATTED')
      READ (1) ( U(I) , I = 1 , NNODE )
      RETURN
      END
C
C
      SUBROUTINE PLTSAI ( DS, NSTEP, START, CRD, S )
      IMPLICIT REAL*8 ( A-H , O-Z )
      DIMENSION CRD(2,3), X(3), Y(3), S(3)
      IF ( NSTEP .EQ. 0 ) RETURN
      DO I = 1 , 3
      X(I) = CRD(1,I)
      Y(I) = CRD(2,I)
      END DO
      DO LEVEL = 1 , NSTEP
      SXY = START + (LEVEL-1) * DS
      K = 1
      DO I = 1 , 3
      J = I + 1
      IF ( I .EQ. 3 ) J = 1
      IF ( (S(I)-SXY)*(S(J)-SXY) .LT. 0 ) THEN
           T=(SXY-S(I))/(S(J)-S(I))
           X0=X(J)*T+(1.D0-T)*X(I)
           Y0=Y(J)*T+(1.D0-T)*Y(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 TEMPA (MXE,MXN,ND,NE,NNODE,NODEX,XCOORD,YCOORD,
     *   S,B,U,NSTEP)
      IMPLICIT REAL*8 ( A-H , O-Z )
      DIMENSION S(ND),NODEX(MXE,ND),XCOORD(MXN),YCOORD(MXN), 
     * B(2,ND),U(MXN)
      CHARACTER FILENAME*12
      FILENAME = "TEMPERAT.OUT"
      CALL PLTLGO ( FILENAME )
      CALL CONTOUR (MXE,MXN,ND,NE,NNODE,NODEX,XCOORD,YCOORD,
     * S,B,U,NSTEP)
      CALL PLTEXT
      RETURN
      END
C
C
      SUBROUTINE CONTOUR (MXE,MXN,ND,NE,NNODE,NODEX,XCOORD,YCOORD,
     *   S,B,P,NSTEP)
      IMPLICIT REAL*8 ( A-H , O-Z )
      DIMENSION S(ND),NODEX(MXE,ND),XCOORD(MXN),
     *   YCOORD(MXN),B(2,ND),P(MXN)
C
      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,YCOORD)
      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
      CALL PLTSAI ( DS, NSTEP, PPMIN, B, S )
      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 = DMIN1 ( QMIN , Q(I) )
      QMAX = DMAX1 ( QMAX , Q(I) )
      END DO
      RETURN
      END
C
C
      SUBROUTINE BOUND (MXE,MXN,ND,NE,NODEX,XCOORD,YCOORD)
      IMPLICIT REAL*8 ( A-H , O-Z )
      DIMENSION XCOORD(MXN),YCOORD(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(NODEX(IEL,I)) , YCOORD(NODEX(IEL,I)) )
      CALL XDRAW ( XCOORD(NODEX(IEL,J)) , YCOORD(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