PROGRAM SET27STAX3D
C=======================================================================
C THREE DIMENSIONAL DATA GENERATING PROGRAM FOR 3DSTATIC8QFXCOMBINE.FOR
C-------------------- TORSION PROBLEM APPLICATION ----------------------
C
C BOUNDARY CONDITIONS
C LEFT END: FIXED IN ALL DIRECTIONS
C RIGHT END: TWIST MOMENT IS APPLIED USING FORCES AT 4 POINTS
C ALL OTHERS: NO NORMAL FOECE, NO SHEAR FORCE
C DOMAIN: RECTANGULAR
C ELEMENT: 27-NODED ISOPARAMETRIC HEXA 3D PARABOLIC ELEMENT
C DOMAIN DISCRETIZATION: ALL ELEMENTS ARE THE SAME SIZE.
C EIJI FUKUMORI JUNE 18, 2010
C=======================================================================
IMPLICIT REAL*8 ( A-H , O-Z )
PARAMETER ( ND=27, INTEPT=3, MXE=22000, MXN=40000, MXB=2000 )
PARAMETER (YOUNG=69.D6,POISSON=0.3D0 )
C=======================================================================
DIMENSION NODEX(MXE,ND),
* XCOORD(MXN), YCOORD(MXN) ,ZCOORD(MXN),
* IBNDFX(MXB), IBNDFY(MXB) ,IBNDFZ(MXB),
* BVX(MXB) , BVY(MXB) ,BVZ(MXB) ,
* IBFORCEX(MXB), IBFORCEY(MXB), IBFORCEZ(MXB),
* BVFORCEX(MXB), BVFORCEY(MXB), BVFORCEZ(MXB)
DIMENSION NEUTRAL(MXN), NODECTOP(MXN)
CHARACTER INPFILE*14,EXFILE*3
LOGICAL YES
C=======================================================================
DATA INPFILE / 'STATIC3D27.DAT' /
C=======================================================================
F(X) = X
G(X) = X
C=======================================================================
TLX = 20.D0
TLY = 2.D0
TLZ = 2.D0
NEX = 20
NEY = 5
NEZ = 5
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=======================================================================
WRITE (*,*)'YOUNG MODULUS = ',YOUNG
WRITE (*,*)'POISSON RATIO = ',POISSON
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
CALL BANDWID ( ND, MXE, NE, NODEX , NBW )
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=======================================================================
C BOUNDARY CONDITIONS
C==== FIRST FIND
C--------- NAVIER EQUATIONS
NBFX = 0
NBFY = 0
NBFZ = 0
NFORCEX = 0
NFORCEY = 0
NFORCEZ = 0
C-----------------------------------------------
C--------- X DIRECTION ON FACE -X ALL NODES
DO K = 1 , NDZ
DO J = 1 , NDY
NBFX = NBFX + 1
IBNDFX(NBFX) = J + (K-1)*NDY
BVX(NBFX) = 0.D0
END DO
END DO
C-----------------------------------------------
C--------- Y DIRECTION ON FACE -X ALL NODES
DO K = 1 , NDZ
DO J = 1 , NDY
NBFY = NBFY + 1
IBNDFY(NBFY) = J + (K-1)*NDY
BVY(NBFY) = 0.D0
END DO
END DO
C-----------------------------------------------
C--------- Z DIRECTION ON FACE -X ALL NODES
DO K = 1 , NDZ
DO J = 1 , NDY
NBFZ = NBFZ + 1
IBNDFZ(NBFZ) = J + (K-1)*NDY
BVZ(NBFZ) = 0.D0
END DO
END DO
C-----------------------------------------------
WRITE(*,*) 'NUMBER OF DIRICHLET BOUNDARY CONDITIONS IN X=',NBFX
WRITE(*,*) 'NUMBER OF DIRICHLET BOUNDARY CONDITIONS IN Y=',NBFY
WRITE(*,*) 'NUMBER OF DIRICHLET BOUNDARY CONDITIONS IN Z=',NBFZ
C===============================================
C================ SECOND KIND
NFORCEX = 0
NFORCEY = 0
NFORCEZ = 0
C--------- FACE OF +X
FORCEINZ = 1000.D0
FORCEINY = 1000.D0
C
C--------- RIGHT END FACE
C
C <---F
C ----3---- A
C I I I
C I I I
C F 2 4 F
C I I I
C I I I
C V ----1----
C F--->
C
C
C----------- POINT 1 ---------
NFORCEY = NFORCEY + 1
IBFORCEY(NFORCEY) = NDY*NDZ*(NDX-1) + 1 + NEY/2*2
BVFORCEY(NFORCEY) = FORCEINY
C----------- POINT 3 ---------
NFORCEY = NFORCEY + 1
IBFORCEY(NFORCEY) = NDY*NDZ*(NDX-1) + 1 + NDY*(NDZ-1) + NEY/2*2
BVFORCEY(NFORCEY) = -FORCEINY
C----------- POINT 2 ---------
NFORCEZ = NFORCEZ + 1
IBFORCEZ(NFORCEZ) = NDY*NDZ*(NDX-1) + 1 + NDY*(NEZ/2*2)
BVFORCEZ(NFORCEZ) = -FORCEINZ
C----------- POINT 4 ---------
NFORCEZ = NFORCEZ + 1
IBFORCEZ(NFORCEZ) = NDY*NDZ*(NDX-1) + 1 + NDY*(NEZ/2*2) + NDY-1
BVFORCEZ(NFORCEZ) = FORCEINZ
C-----------------------------------------------
WRITE(*,*) 'NUMBER OF NEWMANN BOUNDARY CONDITIONS IN X=',NFORCEX
WRITE(*,*) 'NUMBER OF NEWMANN BOUNDARY CONDITIONS IN Y=',NFORCEY
WRITE(*,*) 'NUMBER OF NEWMANN BOUNDARY CONDITIONS IN Z=',NFORCEZ
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,*) 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), ZCOORD(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
WRITE(1,*) NBFZ
DO I = 1 , NBFZ
WRITE (1,*) IBNDFZ(I), BVZ(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
WRITE(1,*) NFORCEZ
IF ( NFORCEZ .GT. 0 ) THEN
DO I = 1 , NFORCEZ
WRITE(1,*) IBFORCEZ(I), BVFORCEZ(I)
END DO
ENDIF
C---------- FINAL
CLOSE (1)
C------------------------ NEUTRAL NODAL NUMBERS ------------------------
C NODEL NUMBERS ALONG NEUTRAL LINE OF BEAM
NNEUT = 1
NEUTRAL(NNEUT) = 2*NEZ/2*NDY + 2*NEY/2 + 1
DO I = 1 , NDX-1
NNEUT = NNEUT + 1
NEUTRAL(NNEUT) = NEUTRAL(NNEUT-1) + NDY*NDZ
END DO
C
OPEN ( 1, FILE='NEUTRAL.DAT', STATUS = 'UNKNOWN')
WRITE (1,*) NNEUT
DO I = 1 , NNEUT
WRITE (1,*) NEUTRAL(I)
END DO
CLOSE (1)
C
C------------- NODES ALONG THE CENTER OF THE TOP FACE(+Z) --------------
NCTOP = 1
NODECTOP(NCTOP) = 1 + NDY*(NDZ-1) + NEY/2*2
DO I = 1 , 2*NEX+1
NCTOP = NCTOP + 1
NODECTOP(NCTOP) = NODECTOP(NCTOP-1) + NDY*NDZ
END DO
C
OPEN ( 1, FILE='NODECTOP.DAT', STATUS = 'UNKNOWN')
WRITE (1,*) NCTOP
DO I = 1 , NCTOP
WRITE (1,*) NODECTOP(I)
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=',NBFX+NBFY+NBFZ,' )'
CLOSE (1)
C=======================================================================
OPEN ( 1, FILE='ELEMENT3D27ND.DAT', STATUS='UNKNOWN' )
C
DO IEL = 1 , NEY*NEZ
DO J = 1 , 4
K = 2*J-1
L = 2*J
M = 2*J+1
IF ( J .EQ. 4 ) M = 1
WRITE(1,*) YCOORD(NODEX(IEL,K)), ZCOORD(NODEX(IEL,K))
WRITE(1,*) YCOORD(NODEX(IEL,L)), ZCOORD(NODEX(IEL,L))
WRITE(1,*) YCOORD(NODEX(IEL,M)), ZCOORD(NODEX(IEL,M))
WRITE(1,*)
END DO
C
DO J = 1 , 4
K = 2*J-1+18
L = 2*J + 18
M = 2*J+1+18
IF ( J .EQ. 4 ) K = 19
WRITE(1,*) YCOORD(NODEX(IEL,K)), ZCOORD(NODEX(IEL,K))
WRITE(1,*) YCOORD(NODEX(IEL,L)), ZCOORD(NODEX(IEL,L))
WRITE(1,*) YCOORD(NODEX(IEL,M)), ZCOORD(NODEX(IEL,M))
WRITE(1,*)
END DO
C
DO J = 1 , 4
K = 2*J-1
L = 2*J + 8
M = 2*J+17
WRITE(1,*) YCOORD(NODEX(IEL,K)), ZCOORD(NODEX(IEL,K))
WRITE(1,*) YCOORD(NODEX(IEL,L)), ZCOORD(NODEX(IEL,L))
WRITE(1,*) YCOORD(NODEX(IEL,M)), ZCOORD(NODEX(IEL,M))
WRITE(1,*)
END DO
C
DO J = 9 , 17 , 2
WRITE(1,*) YCOORD(NODEX(IEL,J)), ZCOORD(NODEX(IEL,J))
WRITE(1,*) YCOORD(NODEX(IEL,J)), ZCOORD(NODEX(IEL,J))
WRITE(1,*)
END DO
C
WRITE(1,*) YCOORD(NODEX(IEL,18)), ZCOORD(NODEX(IEL,18))
WRITE(1,*) YCOORD(NODEX(IEL,18)), ZCOORD(NODEX(IEL,18))
WRITE(1,*)
WRITE(1,*) YCOORD(NODEX(IEL,27)), ZCOORD(NODEX(IEL,27))
WRITE(1,*) YCOORD(NODEX(IEL,27)), ZCOORD(NODEX(IEL,27))
WRITE(1,*)
END DO
CLOSE (1)
C=======================================================================
OPEN ( 1, FILE='FORCES3D27ND.DAT', STATUS='UNKNOWN' )
WRITE (1,*) 'FORCE-DIRECTION YCOORD ZCOORD FORCE-STRENGTH'
DO I = 1 , NFORCEY
NODE = IBFORCEY(I)
WRITE (1,*) 'FORCE-Y', YCOORD(NODE), ZCOORD(NODE), BVFORCEY(I)
WRITE (1,*) 'FORCE-Y', YCOORD(NODE), ZCOORD(NODE), BVFORCEY(I)
WRITE (1,*)
END DO
DO I = 1 , NFORCEZ
NODE = IBFORCEZ(I)
WRITE (1,*) 'FORCE-Z', YCOORD(NODE), ZCOORD(NODE), BVFORCEY(I)
WRITE (1,*) 'FORCE-Z', YCOORD(NODE), ZCOORD(NODE), BVFORCEY(I)
WRITE (1,*)
END DO
CLOSE (1)
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