PROGRAM MAKING_TEST_DATA27
C=======================================================================
C                 27-NODED HEXA PARABOLIC ELEMENT
C   THREE DIMENSIONAL DATA GENERATING PROGRAM FOR LANCZOS-SOUNDX.FOR
C    PROJECT NAME: EIGENVALUE
C    BOUNDARY CONDITIONS
C            ALL BOUNDARY: NATURAL
C    DOMAIN: RECTANGULAR
C    ELEMENT: 27-NODED ISOPARAMETRIC HEXA 3D PARABOLIC ELEMENT
C     EIJI FUKUMORI DECEMBER 9, 2010
C=======================================================================
      IMPLICIT REAL*8 ( A-H , O-Z )
      PARAMETER ( ND=27, INTEPT=3 )
      PARAMETER ( MXE=220000, MXN=4000000, MXB=200000 )
C=======================================================================
      DIMENSION NODEX(MXE,ND), XCOORD(MXN), YCOORD(MXN) ,ZCOORD(MXN)
      CHARACTER INPFILE*13,EXFILE*3
      CHARACTER INFNAME*12
      LOGICAL YES, RSTRAINT
C=======================================================================
      DATA INFNAME / 'EIGN3D27.DAT' /
C=======================================================================
C------------------------ SPEED OF SOUND -------------------------------
      WSPD = 340.D0
C-----------------------------------------------------------------------
C TLX=LENGTH, TLY=HEIGHT, TLZ=DEPTH
      TLX = 4.D0*DATAN(1.D0)
      TLY = TLX/2.D0
      TLZ = TLX/3.D0
C-----------------------------------------------------------------------
      NELEMENT = 10
      ALMIN = DMIN1 ( TLX, TLY, TLZ )
      NEX = NELEMENT * (TLX/ALMIN + 0.001D0)
      NEY = NELEMENT * (TLY/ALMIN + 0.001D0)
      NEZ = NELEMENT * (TLZ/ALMIN + 0.001D0)
      WRITE (*,*) NEX, NEY, NEZ
C=======================================================================
      DX=TLX/(2*NEX)
      DY=TLY/(2*NEY)
      DZ=TLZ/(2*NEZ)
      NDX = 2*NEX+1
      NDY = 2*NEY+1  
      NDZ = 2*NEZ+1 
C=======================================================================
C      NEZ: NEMBER OF VERTICAL   ELEMENTS
C      NEY: NEMBER OF DEPTH      ELEMENTS
C      NEX: NEMBER OF HORIZONTAL ELEMENTS
C=======================================================================
C--------- EIGRNVALUE PARAMETER
      MXNEIGEN=MIN0 ( 20, NNODE/10 )
      IF ( MXNEIGEN .LT. 9 ) MXNEIGEN = 9
C-----------------------------------------------------------------------
C--------- DELTA = SHIFT PARAMETER THAT MUST BE > 0.
      AL = DMAX1 ( TLX, TLY, TLZ )
      DELTA = (4.D0*DATAN(1.D0)/AL)**2 * 0.9
C=======================================================================
C                    ELEMENT CREATION
      NE = 0
      DO I = 1 , NEX
      DO J = 1 , NEZ
      DO K = 1 , NEY
      NE = NE + 1
      IF ( NE .GT. MXE ) STOP 'NE > MXE'
C------- BOTTOM FLOOR
      NODEX(NE,1)  = (I-1)*NDY*NDZ*2 + (J-1)*NDY*2 + (K-1)*2 + 1
      NODEX(NE,2)  = NODEX(NE,1) + NDY*NDZ
      NODEX(NE,3)  = NODEX(NE,2) + NDY*NDZ
      NODEX(NE,4)  = NODEX(NE,3) + 1
      NODEX(NE,5)  = NODEX(NE,4) + 1
      NODEX(NE,8)  = NODEX(NE,1) + 1
      NODEX(NE,7)  = NODEX(NE,8) + 1
      NODEX(NE,9)  = NODEX(NE,2) + 1
      NODEX(NE,6)  = NODEX(NE,9) + 1
C------- MIDIUM FLOOR
      NODEX(NE,10) = NODEX(NE,1) + NDY
      NODEX(NE,17) = NODEX(NE,10) + 1
      NODEX(NE,16) = NODEX(NE,17) + 1
      NODEX(NE,11) = NODEX(NE,10) + NDY*NDZ
      NODEX(NE,18) = NODEX(NE,11) + 1
      NODEX(NE,15) = NODEX(NE,18) + 1
      NODEX(NE,12) = NODEX(NE,11) + NDY*NDZ
      NODEX(NE,13) = NODEX(NE,12) + 1
      NODEX(NE,14) = NODEX(NE,13) + 1
C------- TOP FLOOR
      NODEX(NE,19) = NODEX(NE,10) + NDY
      NODEX(NE,26) = NODEX(NE,19) + 1
      NODEX(NE,25) = NODEX(NE,26) + 1
      NODEX(NE,20) = NODEX(NE,19) + NDY*NDZ
      NODEX(NE,27) = NODEX(NE,20) + 1
      NODEX(NE,24) = NODEX(NE,27) + 1
      NODEX(NE,21) = NODEX(NE,20) + NDY*NDZ
      NODEX(NE,22) = NODEX(NE,21) + 1
      NODEX(NE,23) = NODEX(NE,22) + 1
      END DO
      END DO
      END DO
      WRITE (*,*) 'NUMBER OF ELEMENTS =', NE
C=======================================================================
C--------0---------0---------0---------0---------0---------0---------0--
C                 NODAL COORDINATE CREATION
      NNODE = 0
      DO I = 1 , NDX
      DO K = 1 , NDZ
      DO J = 1 , NDY
      NNODE = NNODE + 1
      IF ( NNODE .GT. MXN ) STOP 'NNODE > MXN'
      XCOORD(NNODE) = (I-1)*DX
      YCOORD(NNODE) = (J-1)*DY
      ZCOORD(NNODE) = (K-1)*DZ
      END DO
      END DO
      END DO
      WRITE (*,*) 'NUMBER OF NODAL POINTS =', NNODE
C=======================================================================
      IWW = 2
      OPEN  (IWW, FILE='DOMAIN.SUM', STATUS='UNKNOWN')
      WRITE (IWW,*) 'SPEED OF SOUND =', WSPD
      WRITE (IWW,*) 'TLX, TLY, TLZ =', TLX, TLY, TLZ
      WRITE (IWW,*) 'TOTAL VOLUME =', TOTALVOL
      WRITE (IWW,*) 'NEX, NEY, NEZ =', NEX, NEY, NEZ
      WRITE (IWW,*) 'MAXIMUM SURFACE LENGTH =', AL
      WRITE (IWW,*) 'SHIFT PARAMETER (DELTA) =', DELTA
      WRITE (IWW,*) 'MAX EIGEN (ALPHA SQ) =', EIGENMX
      WRITE (IWW,*) 'NE =', NE
      WRITE (IWW,*) 'NNODE =', NNODE
      WRITE (IWW,*) 'MXNEIGEN =', MXNEIGEN
      CLOSE (IWW)
C=======================================================================
C                     CREATING DATA FILES
      IW = 1
      OPEN (IW, FILE=INFNAME, STATUS='UNKNOWN')
      WRITE (IW,*) DELTA, MXNEIGEN, WSPD
C---------- ELEMENTS
      WRITE(IW,*) NE
      DO I = 1 , NE
      WRITE (IW,*) I, (NODEX(I,J), J = 1 , ND )
      END DO
C---------- COORDINATES OF NONAL POINTS
      WRITE(IW,*) NNODE
      DO I = 1 , NNODE
      WRITE(IW,*) I,XCOORD(I), YCOORD(I), ZCOORD(I)
      END DO
C=======================================================================
C---------- FINAL
      CLOSE (IW)
C=======================================================================
C------ CREATION OF PARAMETER FILE TO BE USED IN INCLUDE STATEMENT
      CALL BANDWID ( ND, MXE, NE, NODEX , NBW  )
      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 ( MXENGN=',MXNEIGEN,')'
      CLOSE (1)
C=======================================================================
      STOP 'NORMAL TERMINATION'
      END
C
C
      SUBROUTINE BANDWID (  ND, MXE, NE, NODEX , NBW  )
      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 (*,*) 'HALH BANDWIDTH =', NBW
      RETURN
      END