PROGRAM ELECHECK
C=======================================================================
C 3-DIM ELEMENT CHECKER
C ELEMENT : 8-NODED ISO-PARAMETERIC
C NODE NUMBERING ORDER: (+1,-1,-1),(+1,+1,-1),(-1+1,-1),(-1,-1,-1),
C (+1,-1,+1),(+1,+1,+1),(-1+1,+1),(-1,-1,+1)
C CHECK ITEMS:
C GLOBAL STIFFNESS MATRIX'S BANDWIDTH
C ELEMENT VOLUE
C DETERMINANT OF JACOBIAN AT INTEGRATION POINTS
C ORIGINAL:1984 EIJI FUKUMORI, REVISED 2008
C=======================================================================
IMPLICIT REAL*8 ( A-H , O-Z )
INCLUDE 'PARAM.DAT'
CCCCCCCCC PARAMETER ( ND=8,MXE=1000000,MXN=1000000,INTEPT=2 )
PARAMETER ( INTEPT3=INTEPT**3 )
DIMENSION NODEX(MXE,ND),XCOORD(MXN),YCOORD(MXN),ZCOORD(MXN),
* G(ND), F(ND),BX(3,ND),
* SAI(INTEPT), W(INTEPT),BPP(3,ND,INTEPT,INTEPT,INTEPT),
* SF(ND,INTEPT,INTEPT,INTEPT),
* VOLUME(MXE), DETJACOB(MXE,INTEPT3)
CHARACTER INPFILE*12, OUTFILE*12, EXFILE*3
CHARACTER CK(MXE)*8
LOGICAL YES
WRITE (*,*)' ELEMENTCHECK.FOR AN ELEMENT CHECKER PROGRAM'
IF ( ND .NE. 8 ) STOP 'SORRY NO SERVICE FOR ND=27'
C=======================================================================
C PRE-CALCULATION FOR ISO-PARAMETRIC ELEMENT
CALL GRULE ( INTEPT, SAI, W )
CALL DERIV ( ND, INTEPT, G, F, SAI, BPP )
CALL SHAPEF( ND, INTEPT, G, SAI, SF )
C=======================================================================
C FILE MANAGEMENT
INPFILE = 'EIGN3D8.DAT'
OUTFILE ='ELCHECK.RST'
C=======================================================================
WRITE(*,*)' INPUT FILE NAME =', INPFILE
CALL INPUT ( ND,MXE,MXN,NE,NNODE,DELTA,MXNEIGEN,
* NODEX, XCOORD,YCOORD, ZCOORD, INPFILE )
C=======================================================================
CALL BANDWID ( MXE, ND, NE,NODEX, NBW )
WRITE(*,*) 'BANDWIDTH =', NBW
C=======================================================================
EXFILE ='NEW'
IW = 1
INQUIRE ( FILE=OUTFILE, EXIST=YES )
IF ( YES ) EXFILE='OLD'
OPEN ( IW, FILE=OUTFILE, STATUS=EXFILE,FORM='FORMATTED' )
WRITE(IW,*) 'ELEMENT TYPE = 8-NODED HEXAHEDRAL(FIXED)'
WRITE(IW,*) 'NUMBER OF ELEMENTS =', NE
WRITE(IW,*) 'NUMBER OF NODAL POINTS =', NNODE
WRITE(IW,*) 'BANDWIDTH =', NBW
C=======================================================================
CALL VCOMP ( INTEPT,ND,ND3,MXE,MXN,NE,BPP,SF,W,BX,NODEX,
* XCOORD,YCOORD,ZCOORD,VOLUME, DETJACOB,INTEPT3,CK,TVOL )
C=======================================================================
WRITE (IW,*) 'TOTAL VOLUME =', TVOL
DV = TVOL/NE
DX = DV**(1.D0/3.D0)
EIGENMX = ( (4.D0*DATAN(1.D0))/(2.D0*DX) ) ** 2
WRITE (IW,*) 'MAXIMUM EIGEN VALUE(ALPHA SQUARE) =', EIGENMX
C=======================================================================
NGTVVOL = 0
DO I = 1 , NE
IF ( VOLUME(I) .GT. 0.D0 ) THEN
WRITE (IW,'(I5,1X,A8,9G12.5)') I,CK(I),VOLUME(I)
ELSE
NGTVVOL = NGTVVOL + 1
WRITE (IW,'(I5,1X,A8,9G12.5)') I,CK(I),VOLUME(I),
* (DETJACOB(I,J),J=1,INTEPT3)
END IF
END DO
CLOSE (IW)
C=======================================================================
IF ( NGTVVOL .EQ. 0 ) THEN
WRITE (*,*) 'NOW YOU HAVE DATA FILE FOR EIGENVALUE SOLVER.'
ELSE
WRITE (*,*) 'THERE EXISTS NEGATIVE VOLUME ELEMENT.'
WRITE (*,*) 'NEED TO CHECK ELEMENTS.'
END IF
STOP 'NORMAL TERMINATION'
END
C
C
SUBROUTINE INPUT ( ND,MXE,MXN,NE,NNODE,DELTA,MXNEIGEN,
* NODEX, XCOORD,YCOORD,ZCOORD, INPFILE )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION NODEX(MXE,ND),XCOORD(MXN),YCOORD(MXN),ZCOORD(MXN)
LOGICAL YES
CHARACTER INPFILE*12
IR = 1
INQUIRE ( FILE=INPFILE, EXIST=YES )
IF ( YES ) OPEN ( IR, FILE=INPFILE, STATUS='OLD' )
IF ( .NOT. YES ) STOP' NO INPUT FILE'
C========> PARAMETER
C--- APPROPRIATE VALUES: DELTA=0.1, MXNEIGEN=40
READ (1,*) DELTA, MXNEIGEN, WSPD
C========> ELEMENT
READ (1,*) NE
C------MEMORY CHECK
IF ( NE .GT. MXE ) THEN
WRITE (*,*) 'NE=',NE
WRITE (*,*) 'NE > MXE'
STOP
END IF
DO I = 1 , NE
READ (1,*) IEL, ( NODEX(IEL,J), J = 1 , ND )
END DO
C========> COORDINATES
READ (1,*) NNODE
C------MEMORY CHECK
IF ( NNODE .GT. MXN ) THEN
WRITE (*,*) 'NNODE=',NNODE
WRITE (*,*) 'NNODE > MXN'
STOP
END IF
DO I = 1 , NNODE
READ (1,*) NODE, XCOORD(NODE) , YCOORD(NODE), ZCOORD(NODE)
END DO
CLOSE (1)
RETURN
END
C
C
SUBROUTINE BANDWID ( MXE, ND, 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
RETURN
END
C
C
SUBROUTINE GRULE ( N , SAI , W )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION SAI(N) , W(N)
IF ( N .LT. 2 ) STOP'N<2'
IF ( N .GT. 6 ) STOP'N>6'
IF ( N .EQ. 1 ) THEN
SAI(1) = 0.D0
W(1) = 2.D0
RETURN
END IF
C
IF ( N .EQ. 2 ) THEN
SAI(1) = DSQRT(3.D0)/3.D0
W(1) = 1.D0
END IF
C
IF ( N .EQ. 3 ) THEN
SAI(1) = DSQRT(3.D0/5.D0)
SAI(2) = 0.D0
W(1) = 5.D0/ 9.D0
W(2) = 8.D0/ 9.D0
END IF
C
IF ( N .EQ. 4 ) THEN
SAI(1) = 0.33998104358485D0
SAI(2) = 0.86113631159405D0
W(1) = 0.65214515486255D0
W(2) = 0.34785484513745D0
END IF
C
IF ( N .EQ. 5 ) THEN
SAI(1) = 0.90617984593866D0
SAI(2) = 0.53846931010568D0
SAI(3) = 0.D0
W(1) = 0.23692688505619D0
W(2) = 0.47862867049937D0
W(3) = 5.12D0 / 9.D0
END IF
C
IF ( N .EQ. 6 ) THEN
SAI(1) = 0.23861918608319D0
SAI(2) = 0.66120938646626D0
SAI(3) = 0.93246951420315D0
W(1) = 0.46791393457269D0
W(2) = 0.36076157304814D0
W(3) = 0.17132449237917D0
END IF
C
NN = N / 2
DO I = 1 , NN
J = N - I + 1
SAI(J) = - SAI(I)
W(J) = W(I)
END DO
RETURN
END
C
C
SUBROUTINE DERIV ( ND, INTEPT, G, F, SAI, BPP )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION SAI(INTEPT),BPP(3,ND,INTEPT,INTEPT,INTEPT),
* F(ND), G(ND)
C
C SAI(ND) AND W(ND) ARE COORDINATES OF INTEGRATION POINTS
C AND THE WEIGHTING FACTORS RESPECTIVELY.
DO K = 1 , INTEPT
E1 = SAI (K)
DO L = 1 , INTEPT
E2 = SAI (L )
DO M = 1 , INTEPT
E3 = SAI(M)
COMPUTATION OF BP(J) = DN(J) / DETA1
CALL ISOPARA ( ND, E1 + 0.5D0 , E2 , E3 , F )
CALL ISOPARA ( ND, E1 - 0.5D0 , E2 , E3 , G )
DO I = 1 , ND
BPP(1,I,K,L,M) = F(I) - G(I)
END DO
COMPUTATION OF BP(J)= DN(J)/DETA2
CALL ISOPARA ( ND, E1 , E2 + 0.5D0 , E3 , F )
CALL ISOPARA ( ND, E1 , E2 - 0.5D0 , E3 , G )
DO I = 1 , ND
BPP(2,I,K,L,M) = F(I) - G(I)
END DO
COMPUTATION OF BP(J) = DN(J) / DETA3
CALL ISOPARA ( ND, E1 , E2 , E3 + 0.5D0 , F )
CALL ISOPARA ( ND, E1 , E2 , E3 - 0.5D0 , G )
DO I = 1 , ND
BPP(3,I,K,L,M) = F(I) - G(I)
END DO
C
END DO
END DO
END DO
RETURN
END
C
C
SUBROUTINE SHAPEF ( ND, INTEPT, F, SAI, SF )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION F(ND),SAI(INTEPT),SF(ND,INTEPT,INTEPT,INTEPT)
C
C SAI(ND) IS GAUSS'INTEGRATION COORDINATE.
DO K = 1 , INTEPT
E1 = SAI (K)
DO L = 1 , INTEPT
E2 = SAI( L )
DO M = 1 , INTEPT
E3 = SAI ( M )
CALL ISOPARA ( ND, E1 , E2 , E3 , F )
DO I = 1 , ND
SF(I,K,L,M) = F(I)
END DO
C
END DO
END DO
END DO
RETURN
END
C
C
SUBROUTINE ISOPARA ( ND, E1 , E2 , E3 , F )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION F(ND)
F(1) = 0.125D0*(1.D0- E1)*(1.D0- E2)*(1.D0- E3)
F(2) = 0.125D0*(1.D0+ E1)*(1.D0- E2)*(1.D0- E3)
F(3) = 0.125D0*(1.D0+ E1)*(1.D0+ E2)*(1.D0- E3)
F(4) = 0.125D0*(1.D0- E1)*(1.D0+ E2)*(1.D0- E3)
F(5) = 0.125D0*(1.D0- E1)*(1.D0- E2)*(1.D0+ E3)
F(6) = 0.125D0*(1.D0+ E1)*(1.D0- E2)*(1.D0+ E3)
F(7) = 0.125D0*(1.D0+ E1)*(1.D0+ E2)*(1.D0+ E3)
F(8) = 0.125D0*(1.D0- E1)*(1.D0+ E2)*(1.D0+ E3)
RETURN
END
C
C
SUBROUTINE VCOMP( INTEPT,ND,ND3,MXE,MXN,NE,BPP,SF,W,BX,NODEX,
* XCOORD,YCOORD,ZCOORD,VOLUME, DETJACOB,INTEPT3,CK, TVOL )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION NODEX(MXE,ND),XCOORD(MXN),YCOORD(MXN),ZCOORD(MXN),
* BX(3,ND),BPP(3,ND,INTEPT,INTEPT,INTEPT),
* SF(ND,INTEPT,INTEPT,INTEPT),W(INTEPT), YAC(3,3), VOLUME(MXE),
* DETJACOB(MXE,INTEPT3)
CHARACTER CK(MXE)*8
C
TVOL = 0.D0
DO IEL = 1 ,NE
C
VOLUME (IEL) = 0.D0
DO I = 1 , INTEPT3
CK(IEL)(I:I) = '-'
END DO
C
ICOUNT = 0
DO K = 1 , INTEPT
DO L = 1 , INTEPT
DO M = 1 , INTEPT
WEIGHT = W(K) * W(L) * W(M)
C SEE SUBROUTINE DERIV FOR THE SHAPE FUNCTION DERIVATIVES.
COMPUTATION OF JACOBIAN MATRIX [YAC].
DO I = 1 , 3
DO J = 1 , 3
YAC(I,J) = 0.D0
END DO
END DO
C
DO J = 1 , ND
NODE = NODEX(IEL,J)
DO I = 1 , 3
YAC(I,1) = YAC(I,1) + BPP(I,J,K,L,M) * XCOORD(NODE)
YAC(I,2) = YAC(I,2) + BPP(I,J,K,L,M) * YCOORD(NODE)
YAC(I,3) = YAC(I,3) + BPP(I,J,K,L,M) * ZCOORD(NODE)
END DO
END DO
C THE DETERMINANT,DETJ.
DETJ = YAC(1,1) * ( YAC(2,2)*YAC(3,3)-YAC(2,3)*YAC(3,2) )
* - YAC(1,2) * ( YAC(2,1)*YAC(3,3)-YAC(2,3)*YAC(3,1) )
* + YAC(1,3) * ( YAC(2,1)*YAC(3,2)-YAC(2,2)*YAC(3,1) )
ICOUNT = ICOUNT + 1
DETJACOB(IEL,ICOUNT) = DETJ
IF ( DETJ .LE. 0.D0 ) CK(IEL)(ICOUNT:ICOUNT) = '*'
C------ SOURCE TERM
BETA = WEIGHT * DETJ
DO I = 1 , ND
DO J = 1 , ND
VOLUME(IEL) = VOLUME(IEL) + SF(I,K,L,M) * SF(J,K,L,M) * BETA
END DO
END DO
C
END DO
END DO
END DO
C
TVOL = TVOL + VOLUME(IEL)
END DO
RETURN
END