PROGRAM SETSTA3
C=======================================================================
C      DATA GENERATING PROGRAM FOR STATICX.FOR
C   PROJECT NAME: BEAM000   LEFT END: FIXED    RIGHT END: Y-FIX  X-FREE
C    DOMAIN: RECTANGULAR
C    ELEMENT: FOUR-NODED ISOPARAMETRIC ELEMENT
C    DOMAIN DISCRETIZATION: UNEVEN ELEMENTS WITH VERTICAL SCAN
C     EIJI FUKUMORI DECEMBER 29, 1993 
C=======================================================================
      IMPLICIT REAL*8 ( A-H , O-Z )
      PARAMETER ( ND=4, MXE=2200, MXN=2300, MXB=2000 )
      PARAMETER (YOUNG=3.D5,POISSON=0.2D0, HEIGHT=30.D0,WIDTH=100.D0,
     * TLX = WIDTH,   TLY = HEIGHT,   NEY = 6, NEX = 20,
     * DX = TLX / NEX, DY = TLY / NEY, NDX=NEX+1,NDY=NEY+1)
C=======================================================================
      DIMENSION NODEX(MXE,ND), XCOORD(MXN), YCOORD(MXN),
     * IBNDFX(MXB),IBNDFY(MXB),BVX(MXB),BVY(MXB),
     * IBFORCEX(MXB), IBFORCEY(MXB), BVFORCEX(MXB), BVFORCEY(MXB)
      CHARACTER INPFILE*12,EXFILE*3
      LOGICAL YES
C=======================================================================
      DATA INPFILE / 'STATIC04.DAT' /
C=======================================================================
C      F(X) = X * X * ( 3.- 2.* X )
C      G(X) = X * ( (0.1-1.)* X + (2.- 0.1) )
      F(X) = X 
      G(X) = X
C=======================================================================
C         NEY: NEMBER OF VERTICAL ELEMENTS (NUMBER OF NODES: NEY+1)
C         NEX: NEMBER OF HORIZONTAL ELEMENTS (NUMBER OF NODES: NEX+1)
C  HEIGHT: HEIGHT OF DRIVEN CAVITY FLOW DOMAIN
C=======================================================================
      WRITE (*,*)' YOUNG MODULUS = ',YOUNG
      WRITE (*,*)' POISSON RATIO = ',POISSON
C=======================================================================
C                    ELEMENT CREATION
      NE = 0
      DO J = 1 , NEX
      DO I = 1 , NEY
      NE = NE + 1
      IF ( NE .GT. MXE ) STOP 'NE > MXE'
      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
      END DO
      END DO
C=======================================================================
C                 NODAL COORDINATE CREATION
      NNODE = 0
      DO I = 1 , NDX
      DO J = 1 , NDY
      NNODE = NNODE + 1
      IF ( NNODE .GT. MXN ) STOP 'NNODE > MXN'
      XCOORD(NNODE) = TLX * F( DX*(I-1)/TLX )
      YCOORD(NNODE) = TLY * G( DY*(J-1)/TLY )
      END DO
      END DO
C=======================================================================
C                    BOUNDARY CONDITIONS
C--------- NAVIER EQUATIONS
      NBFX = 0
      NBFY = 0
      NFORCEX = 0
      NFORCEY = 0
C--------- FACE OF -X
      U0 = 0.
      V0 = 0.
      DO J = 1 , NDY
      NBFX = NBFX + 1
      NBFY = NBFY + 1
      IBNDFX(NBFX) = J
      IBNDFY(NBFY) = J
      BVX(NBFX) = U0
      BVY(NBFY) = V0
      END DO
C--------- FACE OF +X
      FORCEINX =  0.
      FORCEINY =  0.
      IF ( (FORCEINX .NE. 0.) .AND. (FORCEINY .NE. 0.) ) THEN
      DO I = 1 , NDY
      NFORCEX = NFORCEX + 1
      NFORCEY = NFORCEY + 1
      IBFORCEX(NFORCEX) = NEX*NDY + I
      IBFORCEY(NFORCEY) = NEX*NDY + I
      IF ( (I.EQ.1) .OR. (I.EQ.NDY) ) THEN
      BVFORCEX(NFORCEX) = FORCEINX/2.
      BVFORCEY(NFORCEY) = FORCEINY/2.
      ELSE
      BVFORCEX(NFORCEX) = FORCEINX
      BVFORCEY(NFORCEY) = FORCEINY
      END IF
      END DO
      END IF
C--------- FACE OF +Y
      FORCEINY = -10.
      IF  (FORCEINY .NE. 0.) THEN
      DO I = 1 , NDX
      NFORCEY = NFORCEY + 1
      IBFORCEY(NFORCEY) = I*NDY
      IF ( (I.EQ.1) .OR. (I.EQ.NDX) ) THEN
      BVFORCEY(NFORCEY) = FORCEINY/2.
      ELSE
      BVFORCEY(NFORCEY) = FORCEINY
      END IF
      END DO
      END IF
C--------- FACE OF -Y
C--------- FIX OF DISPLACEMENT V AT THE RIGHT END.
      V0 = 0.
      NBFY = NBFY + 1
      IBNDFY(NBFY) = NEX*NDY+1
      BVY(NBFY) = V0
C=======================================================================
C                      DATA FILE INQUIRY
      EXFILE = 'NEW'
      INQUIRE ( FILE = INPFILE, EXIST = YES )
      IF ( YES ) EXFILE='OLD'
C=======================================================================
C                     CREATING DATA FILES
C---------- PARAMETERS
      OPEN ( 1, FILE=INPFILE, STATUS = EXFILE )
      WRITE(1,*) YOUNG , POISSON
C---------- ELEMENTS
      WRITE(1,*) NE
      DO I = 1 , NE
      WRITE (1,115) I, (NODEX(I,J), J = 1 , ND )
      END DO
  115 FORMAT ( 5 I 6 )
C---------- COORDINATES OF NONAL POINTS
      WRITE(1,*) NNODE
      DO I = 1 , NNODE
      WRITE(1,*) I,XCOORD(I), YCOORD(I)
      END DO
C---------- DIRICHLET TYPE BOUNDARY CONDITIONS
      WRITE(1,*) NBFX
      DO I = 1 , NBFX 
      WRITE (1,*) IBNDFX(I), BVX(I)
      END DO
      WRITE(1,*) NBFY
      DO  I = 1 , NBFY 
      WRITE (1,*) IBNDFY(I), BVY(I)
      END DO
C=======================================================================
C---------- NUEMANN TYPE BOUNDARY CONDITIONS
      WRITE(1,*) NFORCEX
      IF ( NFORCEX .GT. 0 ) THEN
      DO I = 1 , NFORCEX
      WRITE(1,*) IBFORCEX(I), BVFORCEX(I)
      END DO
      ENDIF
      WRITE(1,*) NFORCEY
      IF ( NFORCEY .GT. 0 ) THEN
      DO I = 1 , NFORCEY
      WRITE(1,*) IBFORCEY(I), BVFORCEY(I)
      END DO
      ENDIF
C---------- FINAL
      CLOSE (1)
C=======================================================================
      CALL BANDWID ( ND, MXE, NE, NODEX , NBW  )
      INTEPT = 2
C------ CREATION OF PARAMETER FILE TO BE USED IN INCLUDE STATEMENT
      OPEN ( 1, FILE='PARAM.DAT', STATUS='UNKNOWN' )
      WRITE (1,*) '      PARAMETER ( ND=',ND,'    )'
      WRITE (1,*) '      PARAMETER ( INTEPT=',INTEPT,')'
      WRITE (1,*) '      PARAMETER ( MXE=',NE,')'
      WRITE (1,*) '      PARAMETER ( MXN=',NNODE,')'
      WRITE (1,*) '      PARAMETER ( MXW=',NBW,')'
      NB = MAX0 (NBFX, NBFY, NFORCEX, NFORCEY )
      WRITE (1,*) '      PARAMETER ( MXB=',NB,'   )'
      CLOSE (1)
C=======================================================================
      STOP "NORMAL TERMINATION"
      END
C
C
      SUBROUTINE BANDWID ( ND, MXE, NE, NODEX , NBW  )
      IMPLICIT REAL*8 ( A-H , O-Z )
      DIMENSION NODEX(MXE,ND)
      NBW = 0
      DO  I = 1 , NE
      DO  J = 1 , ND-1
      DO  K = J+1 , ND
      NBW = MAX0(NBW,IABS(NODEX(I,J)-NODEX(I,K)))
      END DO
      END DO
      END DO
      NBW = NBW + 1
      WRITE(*,*) 'HALF BANDWIDTH + 1 =', NBW
      RETURN
      END