PROGRAM SETHL4
C========= WAVE AROUND CYLINDER ========
C=======================================================================
C      DATA GENERATING PROGRAM FOR HELM4Q.FOR
C   PROJECT: DRIVEN CAVITY FLOW   PROJECT NAME: CYLINDR
C    DOMAIN: RECTANGULER WITH A CYLINDER OBJECT IN FLOW FIELD
C    ELEMENT: FOUR-NODED ISOPARAMETRIC ELEMENT
C    DOMAIN DISCRETIZATION: UNEVEN ELEMENTS WITH VERTICAL SCAN
C     REYNOLDS NUMBER: SPECIFIED IN PARAMETER AS REYNLD
C     VISCOSITY RATIO [VRATIO]: 2ND VISCOSITY / 1ST VISCOSITY
C     EIJI FUKUMORI DECEMBER 29, 1993 
C------------- DIAMETER OF CYLINDER = 1. -------------------------------
C=======================================================================
      IMPLICIT REAL*8 ( A-H , O-Z )
      PARAMETER ( ND=4, MXE=20000, MXN=30000, MXB=3000 )
      PARAMETER ( WAVEH=1.,ALPHA=10.D0,
     * ISIDE = 5, IFRONT = 10, IBACK = 20, NEX = IFRONT + 8 + IBACK,
     * NEY = ISIDE + 8 + ISIDE, DX = 0.4, DY = 0.4,
     * NDX = NEX + 1, NDY = NEY + 1  )
C======================================================================
C
      DIMENSION NODEX(MXE,ND),XCOORD(MXN),YCOORD(MXN),
     * IBNDF(MXB),BV(MXB), ITYPE(MXB)
      DIMENSION ICR(7,32),ICREL(192),ICRND(192),XC(32),YC(32)
      DIMENSION IREL(MXN,10),INOD(MXN),NEW(MXN),JNT(MXN),
     *   NJNT(MXN),XTEMP(MXN),YTEMP(MXN)
      CHARACTER PROJECT*12,EXFILE*3
      LOGICAL YES
C=======================================================================
      DATA PROJECT / 'HEL4DATA.QQQ' /
C------------ COORDINATES OF BASE CYLINDER (-1, -1) - (+1, +1) ---------
      DATA XC /
     *-.707,-.555,-.383,-.195,0.,0.195,0.383,0.555,0.707,0.831,
     *0.924,0.981,1.,0.981,0.924,0.831,0.707,0.555,0.383,0.195,
     *0.,-.195,-.383,-.555,-.707,-.831,-.924,-.981,-1.00,-.981,
     *-.924,-.831 /
      DATA YC /
     *-.707,-.831,-.924,-.981,-1.,-.981,-.924,-.831,-.707,-.555,
     *-.383,-.195,0.,0.195,0.383,0.555,0.707,0.831,0.924,0.981,
     *1.,0.981,0.924,0.831,0.707,0.555,0.383,0.195,0.,-0.195,
     *-0.383,-0.555 /
C=======================================================================
      F(X) = X
      G(X) = X
C=======================================================================
C    C1: ADDED VISCOSITY COMPUTATION PARAMETER
C    C1=0. FOR NO ADDED VISCOSITY, C1=1. FOR UPWIND (OPTIMUM) ADDED 
C    C1 OPTION IS AVALABLE IN NSEQ8DD.FOR SOLVER
C         NEY: NEMBER OF VERTICAL ELEMENTS (NUMBER OF NODES: NEY+1)
C         NEX: NEMBER OF HORIZONTAL ELEMENTS (NUMBER OF NODES: NEX+1)
C  MAXITE: MAXIMUM NUMBER OF ITERATIONS AT EACH TIME INCREMENT
C ITERATION PROCEDURE IS REQUIRED BECAUSE OF NON-LINEAR CONVECTIVE TERM.
C  ERMAX: MAXIMUM DIFFERENCE BETWEEN GUESSED {U&V} AND COMPUTED {U&V}
C=======================================================================
      NE = 0
      DO 10 J = 1 , NEX
      DO 10 I = 1 , NEY
      NE = NE + 1
      NODEX(NE,1) = NDY*(J-1) + I
      NODEX(NE,2) = NODEX(NE,1) + NDY
      NODEX(NE,3) = NODEX(NE,2) + 1
      NODEX(NE,4) = NODEX(NE,1) + 1
   10 CONTINUE
C
      TLX = NEX * DX
      TLY = NEY * DY
      NNODE = 0
      DO 20 I = 1 , NDX
      DO 20 J = 1 , NDY
      NNODE = NNODE + 1
      XCOORD(NNODE) = TLX*F( DX*(I-1)/TLX )
      YCOORD(NNODE) = TLY*G( DY*(J-1)/TLY )
   20 CONTINUE
C======== CENTER OF CYLINDER ========
      I = NDY*(IFRONT+4) + ISIDE + 4 + 1
      XCENTER = XCOORD (I)
      YCENTER = YCOORD (I)
C======== NODES AND ELEMENTS AROUND CYLINDER =======
      K = 0
      DO 71 I = 1 , 8
      K = K + 1
      ICR(7,K) = NDY*(IFRONT+I-1) + ISIDE + 1
   71 CONTINUE
      DO 72 I = 1 , 8
      K = K + 1
      ICR(7,K) = NDY*(IFRONT+8) + ISIDE + I
   72 CONTINUE
      DO 73 J = 1 , 8
      I = 8 - J + 1
      K = K + 1
      ICR(7,K) = NDY*(IFRONT+I) + ISIDE + 8 + 1
   73 CONTINUE
      DO 74 J = 1 , 8
      I = 8 - J + 1
      K = K + 1
      ICR(7,K) = NDY*IFRONT + ISIDE + 1 + I
   74 CONTINUE
      K = 0
      DO 75 I = 1 , 7
      DO 75 J = 1 , 7
      K = K + 1
      ICRND(K) = NDY*(IFRONT+1+I-1) + ISIDE + 1 + J
   75 CONTINUE
      DO 76 I = K+1 , 192
      NNODE = NNODE + 1
      ICRND(I) = NNODE
   76 CONTINUE
      K = 0
      DO 77 I = 1 , 6
      DO 77 J = 1 , 32
      K = K + 1
      ICR(I,J) = ICRND(K)
   77 CONTINUE
      K = 0
      DO 78 I = 1 , 8
      DO 78 J = 1 , 8
      K = K + 1
      ICREL(K) = NEY*(IFRONT+I-1) + ISIDE + J
   78 CONTINUE
      DO 79 I = K+1 , 192
      NE = NE + 1
      ICREL(I) = NE
   79 CONTINUE
      WRITE (*,*)' NUMBER OF ELEMENTS =',NE
      WRITE (*,*)' NUMBER OF NODES    =',NNODE
      K = 0
      DO 55 I = 1 , 6
      DO 44 J = 1 , (32-1)
      K = K + 1
      IEL = ICREL(K)
      NODEX(IEL,1) = ICR(I+1,J)
      NODEX(IEL,2) = ICR(I+1,J+1)
      NODEX(IEL,3) = ICR(I,J+1)
      NODEX(IEL,4) = ICR(I,J)
   44 CONTINUE
      K = K + 1
      IEL = ICREL(K)
      NODEX(IEL,1) = ICR(I+1,32)
      NODEX(IEL,2) = ICR(I+1,1)
      NODEX(IEL,3) = ICR(I,1)
      NODEX(IEL,4) = ICR(I,32)
   55 CONTINUE
C------- COORDINATES -----
      DO 41 I = 1 , 32
      XCOORD(ICR(1,I)) = 1.25*DX*XC(I) + XCENTER
      YCOORD(ICR(1,I)) = 1.25*DX*YC(I) + YCENTER
   41 CONTINUE
      DO 51 J = 1 , 32
      RX = (XCOORD(ICR(7,J))-XCOORD(ICR(1,J))) / 6.
      RY = (YCOORD(ICR(7,J))-YCOORD(ICR(1,J))) / 6.
      DO 51 I = 2 , 6
      XCOORD(ICR(I,J)) = XCOORD(ICR(1,J)) + RX*(I-1)
      YCOORD(ICR(I,J)) = YCOORD(ICR(1,J)) + RY*(I-1)
   51 CONTINUE
C
C======================================================================
C--------- BOUNDARY CONDITIONS ---------
C............ MOMENTUM EQUATIONS .......
C        FACE OF -X
      NBF = 0
      DO J = 1 , NDY
      NBF = NBF + 1
      IBNDF(NBF) = J
      BV(NBF) = WAVEH
      ITYPE(NBF) = 1
      END DO
C------- FACE OF +X
      V0 = 0.
      DO 29 I = 1 , NDY
      NBF = NBF + 1
      IBNDF(NBF) = NEX*NDY + I
      BV(NBF) = 0.
      ITYPE(NBF) = 1
   29 CONTINUE
C======================================================================
      CALL OPTIM4(NODEX,IREL,INOD,NEW,JNT,NJNT,XCOORD,YCOORD,XTEMP,
     * YTEMP,MXE,MXN,NE,NNODE,6 )
      DO I = 1 , NBF
      IBNDF(I) = NEW(IBNDF(I))
      END DO
C=======================================================================
C                      DATA FILE INQUIRY
      EXFILE = 'NEW'
      INQUIRE ( FILE = PROJECT, EXIST = YES )
      IF ( YES ) EXFILE='OLD'
C=======================================================================
C                     MAKING DATA FILES
      OPEN ( 1, FILE=PROJECT, STATUS = EXFILE )
C------ CONTROL PARAMETERS
      WRITE(1,'(F10.5)') ALPHA
      WRITE(1,'(I5)'  ) NE
C------ ELEMENTS
      DO I = 1 , NE
      WRITE (1,'(5I5)') I, (NODEX(I,J), J = 1 , ND )
      END DO
C------ NODAL COORDINATES
      WRITE(1,'(I5)'  ) NNODE
      WRITE(1,120) ( I,XCOORD(I), YCOORD(I) , I = 1, NNODE )
  120 FORMAT ( I5, 2F10.6 )
C------ BOUNDARY CONDITIONS
      WRITE(1,'(I5)'  ) NBF
      WRITE (1,125) ( IBNDF(I), ITYPE(I),BV(I) , I = 1 , NBF )
  125 FORMAT ( 2I5, F10.5 )
      CLOSE (1)
      STOP
      END
C
C
      SUBROUTINE OPTIM4 ( ICON,IREL, INOD, NEW, JNT, NJNT, X, Y, 
     &  XTEMP, YTEMP, MAXEL, MAXNOD, NEL, NNOD, IW )
      IMPLICIT REAL*8 ( A-H , O-Z )
      DIMENSION ICON(MAXEL,4), IREL(MAXNOD,10), INOD(MAXNOD), 
     & NEW(MAXNOD),JNT(MAXNOD),NJNT(MAXNOD),X(MAXNOD),Y(MAXNOD),
     &  XTEMP(MAXNOD), YTEMP(MAXNOD),JTEMP(10), ITEMP(4)
C 
      CALL SETUP( ICON,IREL,INOD,MAXEL,MAXNOD,NEL,NNOD,IBWOR ) 
      IMP = 0 
      CALL OPTIM(IREL,INOD,JNT,NJNT,NEW,MAXNOD,NNOD,IBWOR,MAXDIF,IMP) 
      IBW = MAXDIF
      DO 10 I = 1, NNOD 
      JEND = INOD(I)
      DO  20 J = 1, JEND
      JR = JEND - J + 1 
      JTEMP(J) = IREL(I,JR)
   20 CONTINUE
      DO 10 K = 1, JEND 
      IREL(I,K) = JTEMP(K) 
   10 CONTINUE
      CALL OPTIM(IREL,INOD,JNT,NJNT,NEW,MAXNOD,NNOD,IBW,MAXDIF,IMP) 
      IF ( IMP .LE. 0 ) GO TO  90 
      IBOLD = IBWOR + 1 
      MAXBWH = MAXDIF + 1 
      WRITE(IW,200) IBOLD, MAXBWH 
      DO  30 I = 1, NEL 
      DO  40 J = 1, 4 
      ITEMP(J) = ICON(I,J)
   40 CONTINUE
      DO  30 J = 1, 4 
      ICON(I,J) = NEW( ITEMP(J) ) 
   30 CONTINUE
      DO  50 I = 1, NNOD
      XTEMP(I) = X(I) 
      YTEMP(I) = Y(I) 
   50 CONTINUE
      DO  60 I = 1, NNOD
      N = NEW(I)
      X(N) = XTEMP(I) 
      Y(N) = YTEMP(I) 
   60 CONTINUE
C     WRITE(IW,300) 
      DO  70 I = 1, NNOD
C     WRITE(IW,400) I, NEW(I) 
   70 CONTINUE
      GO TO 100 
   90 IBOLD = IBWOR + 1 
      WRITE(IW,500) IBOLD 
  100 RETURN
  200 FORMAT ( /,6X,'HALF-BANDWIDTH HAS BEEN REDUCED FROM',I4,' TO ', 
     &  I4,/ )
  300 FORMAT ( /,6X,'NEW OPTIMIZED NODE NUMBERING:',//, 
     &  10X,'OLD NODE NO.',10X,'NEW NODE NO.',/ ) 
  400 FORMAT ( 14X,I4,4X,'----->',4X,I4 ) 
  500 FORMAT ( 6X,'NODE NUMBERING CANNOT BE IMPROVED.',//,
     &  6X,'HALF-BANDWIDTH REMAINS AT ',I4 )
      END 
C 
C 
      SUBROUTINE SETUP ( ICON,IREL,INOD,MAXEL,MAXNOD,NEL,NNOD,IBWOR ) 
      DIMENSION ICON(MAXEL,4), IREL(MAXNOD,10), INOD(MAXNOD), KPERM(4,3)
      DATA KPERM / 2, 3, 4, 1, 3, 4, 1, 2, 4, 1, 2, 3 / 
      IBWOR = 0 
      DO 100 I = 1, NNOD
      KOUNT = 0 
      DO 200 J = 1, NEL 
      DO 300 K = 1, 4 
      NOD1 = ICON(J,K)
      NOD2 = ICON(J,KPERM(K,1)) 
      NOD3 = ICON(J,KPERM(K,2)) 
      NOD4 = ICON(J,KPERM(K,3)) 
      IBW = MAX0 ( IABS(NOD1-NOD2), IABS(NOD1-NOD3), IABS(NOD1-NOD4), 
     &  IABS(NOD2-NOD3), IABS(NOD2-NOD4), IABS(NOD3-NOD4) ) 
      IF ( I .NE. NOD1 ) GO TO 300
      IF ( KOUNT .EQ. 0 ) GO TO 11
      DO 10 L = 1, KOUNT
      IF ( IREL(I,L) .EQ. NOD2 ) GO TO 12 
   10 CONTINUE
   11 KOUNT = KOUNT + 1 
      IREL(I,KOUNT) = NOD2
   12 DO 13 L = 1, KOUNT
      IF ( IREL(I,L) .EQ. NOD3 ) GO TO 14 
   13 CONTINUE
      KOUNT = KOUNT + 1 
      IREL(I,KOUNT) = NOD3
   14 DO 15 L = 1, KOUNT
      IF ( IREL(I,L) .EQ. NOD4 ) GO TO 300
   15 CONTINUE
      KOUNT = KOUNT + 1 
      IREL(I,KOUNT) = NOD4
  300 CONTINUE
      IBWOR = MAX0 ( IBW, IBWOR ) 
  200 CONTINUE
      INOD(I) = KOUNT 
  100 CONTINUE
      RETURN
      END 
C 
C 
      SUBROUTINE OPTIM ( IREL, INOD, JNT, NJNT, NEW, MAXNOD, NNOD, IDIFF
     &  ,MAXDIF, IMP )
      DIMENSION IREL(MAXNOD,10),INOD(MAXNOD),JNT(MAXNOD),NJNT(MAXNOD),
     &  NEW(MAXNOD) 
      MAXDIF = IDIFF
      DO 60 IX = 1, NNOD
      DO 20 J = 1, NNOD 
      JNT(J) = 0
      NJNT(J) = 0 
   20 CONTINUE
      MAX = 0 
      I = 1 
      NJNT(1) = IX
      JNT(IX) = 1 
      K = 1 
   30 JEND = INOD(NJNT(I))
      DO 40 JJ = 1, JEND
      KX = IREL(NJNT(I),JJ) 
      IF ( JNT(KX) .GT. 0 ) GO TO 40
      K = K + 1 
      NJNT(K) = KX
      JNT(KX) = K 
      NDIFF = IABS ( I - K )
      IF ( NDIFF .GE. MAXDIF ) GO TO 60 
      MAX = MAX0 ( NDIFF, MAX ) 
   40 CONTINUE
      IF ( K .EQ. NNOD ) GO TO 50 
      I = I + 1 
      GO TO 30
   50 MAXDIF = MAX
      DO 70 J = 1, NNOD 
      NEW(J) = JNT(J) 
   70 CONTINUE
      IMP = IMP + 1 
   60 CONTINUE
      RETURN
      END