PROGRAM SETFEM4
C=======================================================================
C                  DATA GENERATING PROGRAM FOR FEM4Q.FOR
C                           DOMAIN : TLX BY TLY
C                           BOUNDARY CONDITIONS
C                    ON FACES OF -X AND -Y: T(X,Y)=0.
C                    ON FACE OF +X: DT/DN = 0
C                    ON FACE OF +Y: T(X,Y) = SIN ( 0.25*PI*X )
C                      EIJI FUKUMORI NOVEMBER 28, 1995
C                               FEB. 13, 2013
C=======================================================================
      IMPLICIT REAL*8 ( A-H , O-Z )
      PARAMETER ( ND=4, MXE=3400, MXN=3500, MXB=2300, INTEPT=2 )
      PARAMETER (   TLX = 2.D0,   TLY = 2.D0 )
C=======================================================================
      DIMENSION NODEX(MXE,ND),XCOORD(MXN),YCOORD(MXN),IBNDT(MXB),
     *          BVT(MXB), ITYPET(MXB)
      CHARACTER PROJECT*14,EXFILE*3
      LOGICAL YES
C=======================================================================
      DATA PROJECT / 'FEM04INPUT.DAT' /
C=======================================================================
C         NEY: NEMBER OF VERTICAL ELEMENTS (NUMBER OF NODES: NEY+1)
C         NEX: NEMBER OF HORIZONTAL ELEMENTS (NUMBER OF NODES: NEX+1)
C=======================================================================
      WRITE(*,'(44HTYPE IN NUMBER OF ELEMENTS IN X-COORD (NEX)=, $)')
      READ (*,*) NEX
      WRITE(*,'(44HTYPE IN NUMBER OF ELEMENTS IN Y-COORD (NEY)=, $)')
      READ (*,*) NEY
      IF ( NEX .LT. 1 ) NEX = 2
      IF ( NEY .LT. 1 ) NEY = 2
      WRITE(*,*)' NEX =',NEX, '     NEY =', NEY
C=======================================================================
      DX = TLX / NEX
      DY = TLY / NEY
      NDX=NEX+1
      NDY=NEY+1
       EXX = 1.D0
       EXY = 0.D0
       EYY = 1.D0
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) = NDX*(I-1) + 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'
      NODE = NDX*(I-1) + J
      XCOORD(NODE) =  DX*(J-1)
      YCOORD(NODE) =  DY*(I-1)
      END DO
      END DO
C=======================================================================
      PI = 4.D0 * DATAN (1.D0)
C                    BOUNDARY CONDITIONS
C--------- MOMENTUM EQUATIONS AND HEAT EQUATION
      NBT  = 0
C--------- FACE OF +Y AND -Y
      DO J = 1 , NDX
C--------- -Y
      NBT = NBT + 1
      IBNDT(NBT) = J
      ITYPET(NBT) = 1
      BVT(NBT) = 0.D0
C--------- +Y
      NBT = NBT + 1
      IBNDT(NBT) = NDX*(NDY-1) + J
      ITYPET(NBT) = 1
      BVT(NBT) = DSIN ( 0.25D0*PI*XCOORD(IBNDT(NBT)))
      END DO
C--------- FACE OF +X AND -X
C           -X
      DO I = 2 , NEY
      NBT  = NBT  + 1
      IBNDT (NBT ) = NDX*(I-1) + 1
      ITYPET(NBT) = 1
      BVT(NBT ) = 0.D0
      END DO
C=======================================================================
      WRITE (*,*) ' NUMBER OF ELEMENTS (NE) = ',NE
      WRITE (*,*) ' NUMBER OF NODAL POINTS (NNODE) = ',NNODE
C=======================================================================
C                      DATA FILE INQUIRY
      EXFILE = 'NEW'
      INQUIRE ( FILE = PROJECT, EXIST = YES )
      IF ( YES ) EXFILE='OLD'
C=======================================================================
C                     MAKING DATA FILES
C---------- PROJECT
      IR = 1
      OPEN ( IR, FILE=PROJECT, STATUS = EXFILE )
      WRITE (IR,*) EXX, EXY, EYY
      WRITE (IR,*) NE
      DO IEL = 1 , NE
        WRITE (IR,'(10I5)') IEL,(NODEX(IEL,J),J=1,ND)
      END DO
      WRITE (IR,*)  NNODE
      DO NODE = 1 , NNODE
        WRITE (IR,*) NODE,XCOORD(NODE),YCOORD(NODE)
      END DO
      WRITE (IR,*) NBT
      DO I = 1 ,   NBT
        WRITE (IR,*) IBNDT(I), ITYPET(I), BVT(I)
      END DO
      CLOSE (1)
C=======================================================================
      OPEN ( 1, FILE='ELEMENT4.OUT', STATUS = 'UNKNOWN' )
      DO IEL = 1 , NE
      DO J = 1 , ND
      NODE = NODEX(IEL,J)
      WRITE(1,*) XCOORD(NODE),YCOORD(NODE)
      END DO
      NODE = NODEX(IEL,1)
      WRITE(1,*) XCOORD(NODE),YCOORD(NODE)
      WRITE(1,*)
      END DO
      CLOSE (1)
C=======================================================================
      CALL BANDWID ( ND, MXE, NE, NODEX , NBW  )
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,')'
      WRITE (1,*) '      PARAMETER ( MXB=',NBT,')'
      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