PROGRAM GRAPHIC3D8
C========X=========X=========X=========X=========X=========X=========X==
C                   GRAPHIC DATA GENERATION SOFTWARE 
C                                FOR 
C               NSEQ8DD3D8LLDECOMPOSITION-IMPLICIT.FOR
C
C                           1 DECEMBER 1985
C                            EIJI FUKUMORI
C=======================================================================
      PARAMETER ( NF=9, ND=8, MXE=60000, MXN=40000 )
      IMPLICIT REAL*8 ( A-H, O-Z )
      DIMENSION XCOORD(3,MXN), U(3,MXN), NODEX(MXE,ND)
      DIMENSION XPLOT(3,MXN), UPLOT(3,MXN)
      LOGICAL DECISION(MXN)
      CHARACTER FNAME(NF)*11
C=======================================================================
      RATIO = 0.4D0
      INDEX = 3
      VALUE = 0.5
C=======================================================================
      WRITE (*,*) '======== 3D FLOW VECTOR PLOTTER ========'
      WRITE (*,*) 'RATIO=', RATIO
      WRITE (*,*) 'INDEX=', INDEX
      WRITE (*,*) 'VALUE=', VALUE
C=======================================================================
      WRITE (*,*)' READING IN DATA FILES', MXE
      CALL INPUT ( MXE,MXN,ND,NE,NNODE,NODEX,XCOORD,U,NF,FNAME )
      WRITE (*,*)'NE   =', NE
      WRITE (*,*)'NNODE=', NNODE
C=======================================================================
      CALL COMPFACT ( RATIO, MXN, NNODE, XCOORD, U, FACT )
      CALL PLOTUV (INDEX, VALUE, MXE, MXN,ND,NE,NNODE, U, NODEX,
     * XCOORD, XPLOT, UPLOT, NPOINT, FACT, DECISION )
C=======================================================================
      STOP
      END
C
C
      SUBROUTINE PLOTUV ( INDEX, VALUE, MXE,MXN,ND,NE,NNODE,U,NODEX,
     * XCOORD,XPLOT,UPLOT,NPOINT,FACT,DECISION )
      IMPLICIT REAL*8 ( A-H, O-Z )
      PARAMETER ( NSEG=12, AL=0.3D0, BETA=0.08D0 )
      DIMENSION NODEX(MXE,ND), XCOORD(3,MXN), U(3,MXN)
      DIMENSION XPLOT(3,MXN), UPLOT(3,MXN),ISEG(12),JSEG(12),UP(3),
     * XP(3), XMAX(3), XMIN(3)
      LOGICAL DECISION(MXN)
      DATA ISEG / 1,2,3,4,1,2,3,4,5,6,7,8 /
      DATA JSEG / 2,3,4,1,5,6,7,8,6,7,8,5 /
C---------
      NPOINT = 0
      DO IEL = 1 , NE
      DO J = 1 , NSEG
      K = NODEX(IEL,ISEG(J))
      L = NODEX(IEL,JSEG(J))
      DIFF1 = VALUE - XCOORD(INDEX,K)
      DIFF2 = VALUE - XCOORD(INDEX,L)
      IF ( DIFF1*DIFF2 .LE. 0.D0 ) THEN
      CALL LOCATION ( ND,MXN,K,L,INDEX,VALUE,XCOORD,U,UP,XP )
      IF ( DSQRT(UP(1)**2+UP(2)**2+UP(3)**2) .NE. 0.D0 ) THEN
      NPOINT = NPOINT + 1
      XPLOT(1,NPOINT) = XP(1)
      XPLOT(2,NPOINT) = XP(2)
      XPLOT(3,NPOINT) = XP(3)
      UPLOT(1,NPOINT) = UP(1)
      UPLOT(2,NPOINT) = UP(2)
      UPLOT(3,NPOINT) = UP(3)
      END IF
      END IF
      END DO
      END DO
      WRITE (*,*) 'NUMBER OF VELOCITY VECTORS =', NPOINT
C-------------------------------------------
      IF ( NPOINT .LE. 0 ) STOP 'NPOINT <= 0'
C-------------------------------------------
      IF ( INDEX .EQ. 1 ) THEN
      WRITE (*,*) 'YOU ARE ABOUT TO VIEW FLOW VECTOR IN Y-Z PLANE'
      DO I = 1 , NPOINT
      XPLOT(1,I) = XPLOT(2,I)
      XPLOT(2,I) = XPLOT(3,I)
      UPLOT(1,I) = UPLOT(2,I)
      UPLOT(2,I) = UPLOT(3,I)
      END DO
      END IF
      IF ( INDEX .EQ. 2 ) THEN
      WRITE (*,*) 'YOU ARE ABOUT TO VIEW FLOW VECTOR IN X-Z PLANE'
      DO I = 1 , NPOINT
      XPLOT(2,I) = XPLOT(3,I)
      UPLOT(2,I) = UPLOT(3,I)
      END DO
      END IF
      IF ( INDEX .EQ. 3 ) THEN
      WRITE (*,*) 'YOU ARE ABOUT TO VIEW FLOW VECTOR IN X-Y PLANE'
      END IF
C--------------------------------------------------
      CALL SORTING ( MXN, NPOINT, XPLOT, DECISION )
C--------------------------------------------------
      OPEN ( 1, FILE='CUTUV.OUT', STATUS='UNKNOWN' )
      DO I = 1 , NPOINT
      IF ( DECISION(I) ) THEN
      X0 = XPLOT(1,I)
      Y0 = XPLOT(2,I)
      DX = UPLOT(1,I) * FACT
      DY = UPLOT(2,I) * FACT
      RNX =   BETA * DY
      RNY = - BETA * DX
      WRITE (1,*) X0, Y0
      WRITE (1,*) X0+(1.-AL)*DX, Y0+(1.-AL)*DY
      WRITE (1,*) X0+(1.-AL)*DX+RNX, Y0+(1.-AL)*DY+RNY
      WRITE (1,*) X0+DX, Y0+DY
      WRITE (1,*) X0+(1.-AL)*DX-RNX, Y0+(1.-AL)*DY-RNY
      WRITE (1,*) X0+(1.-AL)*DX, Y0+(1.-AL)*DY
      WRITE (1,*)
      END IF
      END DO
      CLOSE (1)
      RETURN
      END
C
C
      SUBROUTINE SORTING ( MXN, NPOINT, XPLOT, DECISION )
      IMPLICIT REAL*8 ( A-H, O-Z )
      DIMENSION XPLOT(3,MXN)
      LOGICAL DECISION(MXN)
C
      DO I = 1 , NPOINT
      DECISION(I) = .TRUE.
      END DO
C
      DO I = 1 , NPOINT-1
      IF ( DECISION(I) ) THEN

      DO J = I+1 , NPOINT
      IF ( XPLOT(1,I) .EQ. XPLOT(1,J) .AND.
     *     XPLOT(2,I) .EQ. XPLOT(2,J)       ) THEN
      DECISION(J) = .FALSE.
      END IF
      END DO

      END IF

      END DO
      RETURN
      END
C
C
      SUBROUTINE LOCATION ( ND,MXN,K,L,INDEX,VALUE,XCOORD,U,UP,XP )
      IMPLICIT REAL*8 ( A-H, O-Z )
      DIMENSION XCOORD(3,MXN), U(3,MXN),UP(3),XP(3)
C----------------------------------------
      DIFF1 = VALUE - XCOORD(INDEX,K)
      DIFF2 = VALUE - XCOORD(INDEX,L)
C----------------------------------------
      IF ( DIFF1 .EQ. 0.D0 ) THEN
      DO J = 1 , 3
      UP(J) = U(J,K)
      XP(J) = XCOORD(J,K)
      END DO
      RETURN
      END IF
C----------------------------------------
      IF ( DIFF2 .EQ. 0.D0 ) THEN
      DO J = 1 , 3
      UP(J) = U(J,L)
      XP(J) = XCOORD(J,L)
      END DO
      RETURN
      END IF
C---------------------------------------------------------------
      T=(VALUE-XCOORD(INDEX,K))/(XCOORD(INDEX,L)-XCOORD(INDEX,K))
      DO J = 1 , 3
      XP(J) = (1.D0-T)*XCOORD(J,K)+T*XCOORD(J,L)
      UP(J) = (1.D0-T)*U(J,K)+T*U(J,L)
      END DO
      RETURN
      END
C
C
C========X=========X=========X=========X=========X=========X=========X==
      SUBROUTINE INPUT ( MXE,MXN,ND,NE,NNODE,NODEX,XCOORD,U,NF,FNAME )
      IMPLICIT REAL*8 ( A-H , O-Z )
      DIMENSION NODEX(MXE,ND),XCOORD(3,MXN),U(3,MXN)
      CHARACTER FNAME(NF)*11, PROJECT*7
      LOGICAL YES
C------- RETURN VALUES: NEARLY ALL VARIBLES IN THIS SUBROUTINE
C========> FILENAME NSDATA.FLN
      INQUIRE ( FILE='NSDATA.FLN', EXIST=YES )
      IF ( .NOT. YES ) STOP' NSDATA.FLN DOES NOT EXIST'
      OPEN ( 1, FILE='NSDATA.FLN', STATUS='OLD' )
      READ(1,*)
      READ(1,'(A7)') PROJECT
      DO I = 1 , NF
      READ(1,'(A11)') FNAME(I)
      END DO
      CLOSE (1)
      WRITE(*,*)' PROJECT NAME:============>>>',PROJECT
C========> FILENAME XXXXXXX.ELE
      OPEN ( 1, FILE=FNAME(2), STATUS='OLD' )
      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
      CLOSE (1)
C========> FILENAME XXXXXXX.NOD
      OPEN ( 1, FILE=FNAME(3), STATUS='OLD' )
      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(J,NODE),J=1,3)
      END DO
      CLOSE (1)
C========> FILENAME XXXXXXX.BIN
      OPEN ( 1, FILE=FNAME(6), STATUS='OLD',FORM='UNFORMATTED' )
      READ (1) TIME
      READ (1) ( U(1,I) , I = 1 , NNODE )
      READ (1) ( U(2,I) , I = 1 , NNODE )
      READ (1) ( U(3,I) , I = 1 , NNODE )
      CLOSE (1)
      RETURN
      END
C
C
      SUBROUTINE COMPFACT ( RATIO,MXN,NNODE,XCOORD,U,FACT )
      IMPLICIT REAL*8 ( A-H , O-Z )
      DIMENSION XCOORD(3,MXN),U(3,MXN),XMAX(3), XMIN(3), DX(3)
C-------------------------------------------------------------
      DO J = 1 , 3
      XMAX(J) = XCOORD(J,1)
      XMIN(J) = XCOORD(J,1)
      END DO
      DO I = 1 , NNODE
      DO J = 1 , 3
      XMAX(J) = DMAX1 ( XMAX(J), XCOORD(J,I) )
      XMIN(J) = DMIN1 ( XMIN(J), XCOORD(J,I) )
      END DO
      END DO
      DO J = 1 , 3
      DX(J) = XMAX(J) - XMIN(J)
      END DO
C----------------------------------------------
      RMAX = DSQRT ( U(1,1)**2+U(2,1)**2+U(3,1)**2 )
      DO I = 1 , NNODE
      RMAX = DMAX1 ( RMAX, DSQRT(U(1,I)**2+U(2,I)**2+U(3,I)**2) )
      END DO
      IF ( RMAX .EQ. 0 ) THEN
      WRITE (*,*) 'NO VELOCITY FIELD EXISTS'
      STOP
      END IF
C---------------------------------------------
      FACT = RATIO * DMAX1 (DX(1), DX(2), DX(3) ) / RMAX
      RETURN
      END