PROGRAM SETAXISYMM4
C=======================================================================
C                 ********** 4-NODED ELEMENT **********
C------------------- ELEMENT NUMBERING: HORIZONTAL SCAN ----------------
C     DATA GENERATING PROGRAM FOR STATIC8QFXAXISYMM-NEW.FOR
C     DOMAIN SIZE: TLX BY TLY
C     BOUNDARY CONDITIONS
C------------ NORMAL DIRECTION
C        ON FACE OF -X: P=P0
C        ON FACE OF +X: P=0
C        ON FACE OF -Y: V=0
C        ON FACE OF +Y: V=0
C------------ TANGENTIAL DIRECTION
C        ON FACE OF -X: NATURAL BOUNDARY
C        ON FACE OF +X: NATURAL BOUNDARY
C        ON FACE OF -Y: NATURAL BOUNDARY
C        ON FACE OF +Y: NATURAL BOUNDARY
C     NOV. 28, 2012
C     EIJI FUKUMORI
C=======================================================================
      IMPLICIT REAL*8 ( A-H , O-Z )
      PARAMETER ( ND=4, INTEPT=2, MXE=31300, MXN=34410, MXB=21000 )
      PARAMETER ( TLX0=0.1D0,TLX1=0.15D0,TLY=0.2D0 )
      PARAMETER ( YOUNG=210.D9, POISSON=0.3D0 )
      PARAMETER ( NEX=11,NEY=22 )
      PARAMETER ( P0=10.D6, ALPHA = 1.D0 )
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)
      DIMENSION NEUTRAL(MXN)
      CHARACTER PROJECT*12,EXFILE*3
      LOGICAL YES
C=======================================================================
      DATA PROJECT / 'STATIC4.DAT' /
C=======================================================================
C         NEX: NEMBER OF ELEMENTS IN RADIAL DIRECTION
C         NEY: NEMBER OF ELEMENTS IN VERTICAL DIRECTION
C         P0, P1, P2: APPLIED HYDRAULIC PRESSURE
C=======================================================================
      THICKNSS = TLX1-TLX0
      DX = THICKNSS / NEX
      DY = TLY / NEY
      NDX = NEX + 1
      NDY = NEY + 1
C=======================================================================
      WRITE (*,*)' YOUNG MODULUS = ',YOUNG
      WRITE (*,*)' POISSON RATIO = ',POISSON
C=======================================================================
C                    ELEMENT CREATION
      NE = 0
      DO I = 1 , NEY
      DO J = 1 , NEX
      NE = NE + 1
      IF ( NE .GT. MXE ) STOP 'NE > MXE'
      NODEX(NE,1) = (I-1)*NDX + J
      NODEX(NE,2) = NODEX(NE,1) + 1
      NODEX(NE,3) = NODEX(NE,2) + NDX
      NODEX(NE,4) = NODEX(NE,1) + NDX
      END DO
      END DO
C=======================================================================
C                 NODAL COORDINATE CREATION
      NNODE = 0
      DO I = 1 , NDY
      DO J = 1 , NDX
      NNODE = NNODE + 1
      IF ( NNODE .GT. MXN ) STOP 'NNODE > MXN'
      XCOORD(NNODE) =  THICKNSS * F(DX*(J-1)/THICKNSS,ALPHA) + TLX0
      YCOORD(NNODE) =  DY*(I-1)
      END DO
      END DO
C=======================================================================
      PI = 4.D0* ATAN (1.D0)
C                    BOUNDARY CONDITIONS
C==== FIRST KIND
C--------- NAVIER EQUATIONS
      NBFX = 0
      NBFY = 0
      NFORCEX = 0
      NFORCEY = 0
C
C------- BOTTOM
      DO I = 1 , NDX
      NBFY = NBFY + 1
      IBNDFY(NBFY) = I
      BVY(NBFY) = 0.D0
      END DO
C------- TOP
      DO I = 1 , NDX
      NBFY = NBFY + 1
      IBNDFY(NBFY) = I + NEY*NDX
      BVY(NBFY) = 0.D0
      END DO
C
C==== SECOND KIND
C--------- FACE OF INNER SIDE
      DF0 = 2.D0*PI*TLX0*DY*P0
C------ RESET MEMORY FOR NEUMANN TYPE BOUNDARY VALUE
      DO I = 1 , NDY
      BVFORCEX(I) = 0.D0
      END DO
C--------- LOAD DISTRIBUTION FOR PARABOLIC ELEMENT UNDER UNIFORM LOAD
      WEIGHT = 1.D0/2.D0
      DO I = 1 , NEY
      NODE = (I-1)*NDX + 1
      NFORCEX = (I-1) + 1
      IBFORCEX(NFORCEX) = NODE
      BVFORCEX(NFORCEX) = BVFORCEX(NFORCEX) + WEIGHT*DF0
      NODE = NODE + NDX
      NFORCEX = I + 1
      IBFORCEX(NFORCEX) = NODE
      BVFORCEX(NFORCEX) =  BVFORCEX(NFORCEX) +WEIGHT*DF0
      END DO
C=======================================================================
      WRITE (*,*) ' NUMBER OF ELEMENTS (NE) = ',NE
      WRITE (*,*) ' NUMBER OF NODAL POINTS (NNODE) = ',NNODE
      WRITE (*,*) ' NUMBER OF DIRICHLET R (NBFX) = ',NBFX
      WRITE (*,*) ' NUMBER OF DIRICHLET Z (NBFY) = ',NBFY
      WRITE (*,*) ' NUMBER OF NEUMANN R (NFORCEX) = ',NFORCEX
      WRITE (*,*) ' NUMBER OF NEUMANN Z (NFORCEY) = ',NFORCEY
C=======================================================================
C                      DATA FILE INQUIRY
      EXFILE = 'NEW'
      INQUIRE ( FILE = PROJECT, EXIST = YES )
      IF ( YES ) EXFILE='OLD'
C=======================================================================
C                     MAKING DATA FILES
C---------- 'PROJECT'.JNK
      IR = 1
C---------- PARAMETERS
      OPEN ( IR, FILE=PROJECT, STATUS = EXFILE )
      WRITE(1,*) YOUNG , POISSON
C---------- ELEMENTS
C---------- ELEMENTS
      WRITE(1,*) NE
      DO I = 1 , NE
      WRITE (1,*) I, (NODEX(I,J), J = 1 , ND )
      END DO
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---------- 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---------- ELEMENT DRAWING
      OPEN ( 1, FILE='ELEMENT4.DAT', STATUS = 'UNKNOWN')
      DO I = 1, NE
      DO J = 1, ND
      WRITE (1,*) XCOORD(NODEX(I,J)), YCOORD(NODEX(I,J))
      END DO
      WRITE (1,*) XCOORD(NODEX(I,1)), YCOORD(NODEX(I,1))
      WRITE (1,*)
      END DO
      CLOSE (1)
C=======================================================================
C------ CREATION OF PARAMETER FILE TO BE USED IN INCLUDE STATEMENT
      CALL BANDWID ( MXE, ND, NE, NODEX, NBW )
      OPEN ( 1, FILE='PARAM.DAT', STATUS='UNKNOWN' )
      WRITE (1,*) '      PARAMETER ( ND=',ND,', INTEPT=',INTEPT,' )'
      WRITE (1,*) '      PARAMETER ( MXE=',NE,', MXN=',NNODE,
     * ', MXB=',MAX0(NBFX,NBFY,NFORCEX,NFORCEY), ', MXW=',NBW, ' )'
      CLOSE (1)
      STOP 'NORMAL TERMINATION'
      END
C
C
      FUNCTION F (X, A)
      IMPLICIT REAL*8 ( A-H , O-Z )
      F = X**A
      RETURN
      END
C
C
      SUBROUTINE BANDWID ( MXE , ND , NE , NODEX , NBW )
      DIMENSION NODEX(MXE,ND)
C------- RETURN VALUE: NBW
      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 (*,*) ' HALH BANDWIDTH =', NBW
      RETURN
      END