PROGRAM GRA12FILE
C=======================================================================
C FINITE ELEMENT GRAPHICS FOR TWO-DIMENSIONAL FEM9 SOLUTIONS
C ELEMENT TYPE: 12-NODED ISOPARAMETRIC ELEMENT
C CREATE FILES CONTAINING XY-COORDINTE VALUES
C ORIGINAL CODE: NSEQ8GRA.FOR, EIJI FUKUMORI, JULY 1985
C 15 FEB. 2013
C=======================================================================
INCLUDE 'PARAM.DAT'
IMPLICIT REAL*8 ( A-H , O-Z )
C ARRAYS
CHARACTER*14 INPFILE
DIMENSION NODEX(MXE,ND), XCOORD(2,MXN), U(MXN),
* F13(ND), F14(ND), F15(ND), F16(ND)
C=======================================================================
WRITE(*,*)' FEM12 GRAPHICS PROGRAM'
WRITE (*,*)' READING IN DATA FILES'
CALL INPUT ( INPFILE,ND,MXE,MXN,NE,NNODE,NODEX,XCOORD,U )
WRITE (*,*)'PROJECT NAME =======>', INPFILE
C=======================================================================
NSTEP = 20
CALL TEMPA (MXE,MXN,ND,NE,NNODE,NODEX,XCOORD,U,NSTEP,
* F13, F14, F15, F16 )
C=======================================================================
STOP 'TERMINATION'
END
C
C
SUBROUTINE INPUT ( INPFILE,ND,MXE,MXN,NE,NNODE,NODEX,XCOORD,U )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION NODEX(MXE,ND),XCOORD(2,MXN), U(MXN)
CHARACTER INPFILE*14
IF ( ND .LE. 2 ) STOP 'NOT ACCEPTABLE ND'
IF ( ND .EQ. 3 ) INPFILE = 'FEM03INPUT.DAT'
IF ( ND .EQ. 4 ) INPFILE = 'FEM04INPUT.DAT'
IF ( ND .EQ. 8 ) INPFILE = 'FEM08INPUT.DAT'
IF ( ND .EQ. 9 ) INPFILE = 'FEM09INPUT.DAT'
IF ( ND .EQ. 12 ) INPFILE = 'FEM12INPUT.DAT'
OPEN ( 1, FILE = INPFILE, STATUS = 'UNKNOWN' )
READ (1,*) EXX, EXY, EYY
READ (1,*) NE
IF ( NE .GT. MXE ) STOP 'NE > MXE'
READ (1,*) (IEL,(NODEX(IEL,J),J=1,ND), I=1,NE)
READ (1,*) NNODE
IF ( NNODE .GT. MXN ) STOP 'NNODE > MXN'
READ (1,*) (NODE,XCOORD(1,NODE),XCOORD(2,NODE),J=1,NNODE)
CLOSE (1)
C=======================================================================
C========> FILENAME XXXXXXX.BIN
OPEN (1,FILE="SOLUTION.BIN",STATUS='UNKNOWN',FORM='UNFORMATTED')
READ (1) ( U(I) , I = 1 , NNODE )
CLOSE (1)
RETURN
END
C
C
SUBROUTINE BOUND (MXE,MXN,ND,NE,NODEX,XCOORD)
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION XCOORD(2,MXN),NODEX(MXE,ND)
LOGICAL LINE
C--------- FOR USE OF ND=12 ELEMENTS ---------------
NDD = 11
DO IEL = 1 , NE
C
DO I = 1 , NDD, 3
LINE = .TRUE.
J = I + 1
K = J + 1
L = K + 1
IF ( I .EQ. 10 ) L = 1
C
DO JEL = 1 , NE
IF ( IEL .NE. JEL ) THEN
DO II = 1 , NDD, 3
JJ = II + 1
KK = JJ + 1
LL = KK + 1
IF ( II .EQ. 10 ) LL = 1
IF ( NODEX(IEL,I).EQ.NODEX(JEL,LL) .AND.
* NODEX(IEL,L).EQ.NODEX(JEL,II) ) THEN
LINE = .FALSE.
EXIT
END IF
END DO
END IF
IF ( .NOT. LINE ) EXIT
END DO
C
IF ( LINE ) THEN
WRITE (1,*) XCOORD(1,NODEX(IEL,I)) , XCOORD(2,NODEX(IEL,I))
WRITE (1,*) XCOORD(1,NODEX(IEL,J)) , XCOORD(2,NODEX(IEL,J))
WRITE (1,*) XCOORD(1,NODEX(IEL,K)) , XCOORD(2,NODEX(IEL,K))
WRITE (1,*) XCOORD(1,NODEX(IEL,L)) , XCOORD(2,NODEX(IEL,L))
WRITE (1,*)
END IF
END DO
C
END DO
RETURN
END
C
C
SUBROUTINE TEMPA (MXE,MXN,ND,NE,NNODE,NODEX,XCOORD,P,NSTEP,
* F13, F14, F15, F16 )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION NODEX(MXE,ND),XCOORD(2,MXN),P(MXN),
* F13(ND), F14(ND), F15(ND), F16(ND)
CHARACTER FILENAME*12
FILENAME = "TEMPEAT.OUT"
CALL PLTLGO ( FILENAME )
CALL CONTOUR ( MXE, MXN,ND,NE,NNODE,NODEX,XCOORD,P,NSTEP,
* F13, F14, F15, F16 )
CALL PLTEXT
RETURN
END
C
C
SUBROUTINE CONTOUR (MXE,MXN,ND,NE,NNODE,NODEX,XCOORD,P,NSTEP,
* F13, F14, F15, F16 )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION NODEX(MXE,ND),XCOORD(2,MXN),P(MXN),S(4),B(2,4),
* IFOUR(4,4), F13(ND), F14(ND), F15(ND), F16(ND)
DATA (IFOUR(1,J),J=1,4) /1, 2, 9, 8 /
DATA (IFOUR(2,J),J=1,4) /2, 3, 4, 9 /
DATA (IFOUR(3,J),J=1,4) /8, 9, 6, 7 /
DATA (IFOUR(4,J),J=1,4) /9, 4, 5, 6 /
C
CALL MINMAX ( MXN, NNODE, P, PPMIN, PPMAX )
IF ( PPMAX .EQ. PPMIN ) RETURN
NSTEP = NSTEP/2*2 + 1
DS = ( PPMAX - PPMIN ) / NSTEP
I = PPMIN/DS
PPMIN= I*DS
C
CALL BOUND ( MXE,MXN,ND,NE,NODEX,XCOORD )
CALL ISOPARA ( ND, -1.D0/3.D0 , -1.D0/3.D0 , F13 )
CALL ISOPARA ( ND, +1.D0/3.D0 , -1.D0/3.D0 , F14 )
CALL ISOPARA ( ND, +1.D0/3.D0 , +1.D0/3.D0 , F15 )
CALL ISOPARA ( ND, -1.D0/3.D0 , +1.D0/3.D0 , F16 )
DO IEL = 1 , NE
X13 = 0.D0
Y13 = 0.D0
P13 = 0.D0
DO J = 1 , ND
X13 = X13 + F13(J)*XCOORD(1,NODEX(IEL,J))
Y13 = Y13 + F13(J)*XCOORD(2,NODEX(IEL,J))
P13 = P13 + F13(J)*P(NODEX(IEL,J))
END DO
C
X14 = 0.D0
Y14 = 0.D0
P14 = 0.D0
DO J = 1 , ND
X14 = X14 + F14(J)*XCOORD(1,NODEX(IEL,J))
Y14 = Y14 + F14(J)*XCOORD(2,NODEX(IEL,J))
P14 = P14 + F14(J)*P(NODEX(IEL,J))
END DO
C
X15 = 0.D0
Y15 = 0.D0
P15 = 0.D0
DO J = 1 , ND
X15 = X15 + F15(J)*XCOORD(1,NODEX(IEL,J))
Y15 = Y15 + F15(J)*XCOORD(2,NODEX(IEL,J))
P15 = P15 + F15(J)*P(NODEX(IEL,J))
END DO
C
X16 = 0.D0
Y16 = 0.D0
P16 = 0.D0
DO J = 1 , ND
X16 = X16 + F16(J)*XCOORD(1,NODEX(IEL,J))
Y16 = Y16 + F16(J)*XCOORD(2,NODEX(IEL,J))
P16 = P16 + F16(J)*P(NODEX(IEL,J))
END DO
C********* 12-NODED ELEMENT DIVIDED INTO 9 SUB-ELEMENTS ************
C------ 1ST SUB-ELEMENT
NODE = NODEX(IEL,1)
B(1,1) = XCOORD(1,NODE)
B(2,1) = XCOORD(2,NODE)
S( 1) = P(NODE)
NODE = NODEX(IEL,2)
B(1,2) = XCOORD(1,NODE)
B(2,2) = XCOORD(2,NODE)
S( 2) = P(NODE)
B(1,3) = X13
B(2,3) = Y13
S( 3) = P13
NODE = NODEX(IEL,12)
B(1,4) = XCOORD(1,NODE)
B(2,4) = XCOORD(2,NODE)
S( 4) = P(NODE)
CALL PLTSAI ( DS, NSTEP, PPMIN, B, S )
C------ 2ND SUB-ELEMENT
NODE = NODEX(IEL,2)
B(1,1) = XCOORD(1,NODE)
B(2,1) = XCOORD(2,NODE)
S( 1) = P(NODE)
NODE = NODEX(IEL,3)
B(1,2) = XCOORD(1,NODE)
B(2,2) = XCOORD(2,NODE)
S( 2) = P(NODE)
B(1,3) = X14
B(2,3) = Y14
S( 3) = P14
B(1,4) = X13
B(2,4) = Y13
S( 4) = P13
CALL PLTSAI ( DS, NSTEP, PPMIN, B, S )
C------ 3RD SUB-ELEMENT
NODE = NODEX(IEL,3)
B(1,1) = XCOORD(1,NODE)
B(2,1) = XCOORD(2,NODE)
S( 1) = P(NODE)
NODE = NODEX(IEL,4)
B(1,2) = XCOORD(1,NODE)
B(2,2) = XCOORD(2,NODE)
S( 2) = P(NODE)
NODE = NODEX(IEL,5)
B(1,3) = XCOORD(1,NODE)
B(2,3) = XCOORD(2,NODE)
S( 3) = P(NODE)
B(1,4) = X14
B(2,4) = Y14
S( 4) = P14
CALL PLTSAI ( DS, NSTEP, PPMIN, B, S )
C------ 4TH SUB-ELEMENT
NODE = NODEX(IEL,12)
B(1,1) = XCOORD(1,NODE)
B(2,1) = XCOORD(2,NODE)
S( 1) = P(NODE)
B(1,2) = X13
B(2,2) = Y13
S( 2) = P13
B(1,3) = X16
B(2,3) = Y16
S( 3) = P16
NODE = NODEX(IEL,11)
B(1,4) = XCOORD(1,NODE)
B(2,4) = XCOORD(2,NODE)
S( 4) = P(NODE)
CALL PLTSAI ( DS, NSTEP, PPMIN, B, S )
C------ 5TH SUB-ELEMENT
K = 1
B(1,K) = X13
B(2,K) = Y13
S( K) = P13
K = 2
B(1,K) = X14
B(2,K) = Y14
S( K) = P14
K = 3
B(1,K) = X15
B(2,K) = Y15
S( K) = P15
K = 4
B(1,K) = X16
B(2,K) = Y16
S( K) = P16
CALL PLTSAI ( DS, NSTEP, PPMIN, B, S )
C------ 6TH SUB-ELEMENT
B(1,1) = X14
B(2,1) = Y14
S( 1) = P14
NODE = NODEX(IEL,5)
B(1,2) = XCOORD(1,NODE)
B(2,2) = XCOORD(2,NODE)
S( 2) = P(NODE)
NODE = NODEX(IEL,6)
B(1,3) = XCOORD(1,NODE)
B(2,3) = XCOORD(2,NODE)
S( 3) = P(NODE)
B(1,4) = X15
B(2,4) = Y15
S( 4) = P15
CALL PLTSAI ( DS, NSTEP, PPMIN, B, S )
C------ 7TH SUB-ELEMENT
NODE = NODEX(IEL,11)
B(1,1) = XCOORD(1,NODE)
B(2,1) = XCOORD(2,NODE)
S( 1) = P(NODE)
K = 2
B(1,K) = X16
B(2,K) = Y16
S( K) = P16
NODE = NODEX(IEL,9)
B(1,3) = XCOORD(1,NODE)
B(2,3) = XCOORD(2,NODE)
S( 3) = P(NODE)
NODE = NODEX(IEL,10)
B(1,4) = XCOORD(1,NODE)
B(2,4) = XCOORD(2,NODE)
S( 4) = P(NODE)
CALL PLTSAI ( DS, NSTEP, PPMIN, B, S )
C------ 8TH SUB-ELEMENT
K = 1
B(1,K) = X16
B(2,K) = Y16
S( K) = P16
K = 2
B(1,K) = X15
B(2,K) = Y15
S( K) = P15
NODE = NODEX(IEL,8)
B(1,3) = XCOORD(1,NODE)
B(2,3) = XCOORD(2,NODE)
S( 3) = P(NODE)
NODE = NODEX(IEL,9)
B(1,4) = XCOORD(1,NODE)
B(2,4) = XCOORD(2,NODE)
S( 4) = P(NODE)
CALL PLTSAI ( DS, NSTEP, PPMIN, B, S )
C------ 9TH SUB-ELEMENT
K = 1
B(1,K) = X15
B(2,K) = Y15
S( K) = P15
NODE = NODEX(IEL,6)
B(1,2) = XCOORD(1,NODE)
B(2,2) = XCOORD(2,NODE)
S( 2) = P(NODE)
NODE = NODEX(IEL,7)
B(1,3) = XCOORD(1,NODE)
B(2,3) = XCOORD(2,NODE)
S( 3) = P(NODE)
NODE = NODEX(IEL,8)
B(1,4) = XCOORD(1,NODE)
B(2,4) = XCOORD(2,NODE)
S( 4) = P(NODE)
CALL PLTSAI ( DS, NSTEP, PPMIN, B, S )
C********************************************************************
END DO
RETURN
END
C
C
SUBROUTINE PLTSAI ( DS, NSTEP, START, CRD, S )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION CRD(2,4), S(4)
IF ( NSTEP .EQ. 0 ) RETURN
DO LEVEL = 1 , NSTEP
SXY = START + (LEVEL-1) * DS
K = 1
DO I = 1 , 4
J = I + 1
IF ( I .EQ. 4 ) J = 1
IF ( (S(I)-SXY)*(S(J)-SXY) .LT. 0 ) THEN
T=(SXY-S(I))/(S(J)-S(I))
X0=CRD(1,J)*T+(1.D0-T)*CRD(1,I)
Y0=CRD(2,J)*T+(1.D0-T)*CRD(2,I)
IF ( K .GT. 0 ) THEN
CALL XMOVE ( X0, Y0 )
ELSE
CALL XDRAW ( X0, Y0 )
END IF
K = -K
END IF
END DO
END DO
RETURN
END
C
C
SUBROUTINE ISOPARA ( ND , E1 , E2 , F )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION F(ND)
F(1) =0.03125D0*(1.D0-E1)*(1.D0-E2)*(-10.D0+9.D0*(E1*E1+E2*E2))
F(2) =0.28125D0*(1.D0-E2)*(1.D0-E1*E1)*(1.D0-3.D0*E1)
F(3) =0.28125D0*(1.D0-E2)*(1.D0-E1*E1)*(1.D0+3.D0*E1)
F(4) =0.03125D0*(1.D0+E1)*(1.D0-E2)*(-10.D0+9.D0*(E1*E1+E2*E2))
F(5) =0.28125D0*(1.D0+E1)*(1.D0-E2*E2)*(1.D0-3.D0*E2)
F(6) =0.28125D0*(1.D0+E1)*(1.D0-E2*E2)*(1.D0+3.D0*E2)
F(7) =0.03125D0*(1.D0+E1)*(1.D0+E2)*(-10.D0+9.D0*(E1*E1+E2*E2))
F(8) =0.28125D0*(1.D0+E2)*(1.D0-E1*E1)*(1.D0+3.D0*E1)
F(9) =0.28125D0*(1.D0+E2)*(1.D0-E1*E1)*(1.D0-3.D0*E1)
F(10)=0.03125D0*(1.D0-E1)*(1.D0+E2)*(-10.D0+9.D0*(E1*E1+E2*E2))
F(11)=0.28125D0*(1.D0-E1)*(1.D0-E2*E2)*(1.D0+3.D0*E2)
F(12)=0.28125D0*(1.D0-E1)*(1.D0-E2*E2)*(1.D0-3.D0*E2)
RETURN
END
C
C
SUBROUTINE MINMAX ( MXN, NNODE, Q, QMIN, QMAX )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION Q(MXN)
QMIN = Q(1)
QMAX = Q(1)
DO I = 1 , NNODE
QMIN = DMIN1 ( QMIN , Q(I) )
QMAX = DMAX1 ( QMAX , Q(I) )
END DO
RETURN
END
C
C======================== GRAPHICS RUTINES =============================
SUBROUTINE PLTLGO ( FILENAME )
IMPLICIT REAL*8 ( A-H , O-Z )
CHARACTER FILENAME*12
OPEN ( 1, FILE=FILENAME, STATUS='UNKNOWN' )
RETURN
END
C
C
SUBROUTINE PLTEXT
IMPLICIT REAL*8 ( A-H , O-Z )
CLOSE (1)
RETURN
END
C
C
SUBROUTINE JCOLOR ( I )
RETURN
END
C
C
SUBROUTINE XMOVE ( X , Y )
IMPLICIT REAL*8 ( A-H , O-Z )
WRITE(1,*) X, Y
RETURN
END
C
C
SUBROUTINE XDRAW ( X , Y )
IMPLICIT REAL*8 ( A-H , O-Z )
WRITE(1,*) X, Y
WRITE(1,*)
RETURN
END