PROGRAM SETFEM12
C=======================================================================
C ********** 12-NODED ELEMENT **********
C------------------- ELEMENT NUMBERING: HORIZONTAL SCAN ----------------
C DATA GENERATING PROGRAM FOR FEMXQ.FOR
C DOMAIN SIZE: TLX BY TLY
C BOUNDARY CONDITIONS
C------------ NORMAL DIRECTION
C ON FACE OF -X: U=0
C ON FACE OF +X: DU/DN=0
C ON FACE OF -Y: U=0
C ON FACE OF +Y: U=SIN((PI/4)*X)
C FEB. 14, 2013
C EIJI FUKUMORI
C=======================================================================
IMPLICIT REAL*8 ( A-H , O-Z )
PARAMETER ( ND=12, INTEPT=3, MXE=31300, MXN=34410, MXB=21000 )
PARAMETER ( TLX1=2.D0,TLY=2.D0 )
C=======================================================================
DIMENSION NODEX(MXE,ND),XCOORD(MXN),YCOORD(MXN),
* IBNDT(MXB),ITYPET(MXB),BVT(MXB)
CHARACTER PROJECT*14
LOGICAL YES
C=======================================================================
DATA PROJECT / 'FEM12INPUT.DAT' /
C=======================================================================
C NEX: NEMBER OF ELEMENTS IN RADIAL DIRECTION
C NEY: NEMBER OF ELEMENTS IN VERTICAL DIRECTION
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=======================================================================
EXX = 1.D0
EXY = 0.D0
EYY = 1.D0
C=======================================================================
DXEDGE = (TLX1/ NEX) / 3.D0
DXCNTR = TLX1/ NEX
DYEDGE = (TLY / NEY) / 3.D0
DYCNTR = TLY / NEY
NDXEDGE = NEX*3 + 1
NDXCNTR = NEX + 1
NDYEDGE = NEY*3 + 1
NDYCNTR = NEY + 1
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) = (I-1)*(NDXEDGE+NDXCNTR*2) + (J-1)*3 + 1
NODEX(NE, 2) = NODEX(NE,1) + 1
NODEX(NE, 3) = NODEX(NE,2) + 1
NODEX(NE, 4) = NODEX(NE,3) + 1
NODEX(NE, 5) = NODEX(NE,4) + (NEX-J)*3 + J + 1
NODEX(NE, 6) = NODEX(NE,5) + NDXCNTR
NODEX(NE, 7) = NODEX(NE,6) + (NEX-J) + J*3 + 1
NODEX(NE, 8) = NODEX(NE,7) - 1
NODEX(NE, 9) = NODEX(NE,8) - 1
NODEX(NE,10) = NODEX(NE,9) - 1
NODEX(NE,11) = NODEX(NE,6) - 1
NODEX(NE,12) = NODEX(NE,5) - 1
END DO
END DO
C=======================================================================
C NODAL COORDINATE CREATION
NNODE = 0
DO I = 1 , NEY
IF ( NNODE .GT. MXN ) STOP 'NNODE > MXN'
C------- LOWER
DO J = 1 , NDXEDGE
NNODE = NNODE + 1
XCOORD(NNODE) = DXEDGE*(J-1)
YCOORD(NNODE) = DYCNTR*(I-1)
END DO
C------- FIRST MID
DO J = 1 , NDXCNTR
NNODE = NNODE + 1
XCOORD(NNODE) = DXCNTR*(J-1)
YCOORD(NNODE) = DYCNTR*(I-1) + DYEDGE
END DO
C------- SECOND MID
DO J = 1 , NDXCNTR
NNODE = NNODE + 1
XCOORD(NNODE) = DXCNTR*(J-1)
YCOORD(NNODE) = DYCNTR*(I-1) + 2.D0*DYEDGE
END DO
END DO
C
C------- TOP
DO J = 1 , NDXEDGE
NNODE = NNODE + 1
XCOORD(NNODE) = DXEDGE*(J-1)
YCOORD(NNODE) = DYCNTR*NEY
END DO
C=======================================================================
PI = 4.D0* ATAN (1.D0)
C BOUNDARY CONDITIONS
C==== FIRST KIND
NBT = 0
C
C------- BOTTOM (-Y)
DO I = 1 , NDXEDGE
NBT = NBT + 1
IBNDT (NBT) = I
ITYPET(NBT) = 1
BVT (NBT) = 0.D0
END DO
C------- TOP (+Y)
DO I = 1 , NDXEDGE
NBT = NBT + 1
IBNDT (NBT) = NNODE - NDXEDGE + I
ITYPET(NBT) = 1
BVT (NBT) = DSIN ( 0.25D0*PI*XCOORD(IBNDT(NBT)))
END DO
C--------- FACE OF LEFT-SIDE (-X)
DO I = 1 , NEY
NODE1 = (I-1)*(NEX*3+1+(NEX+1)*2) + 1
NODE2 = NODE1 + NEX*3 + 1
NODE3 = NODE2 + NEX + 1
NODE4 = NODE3 + NEX + 1
C
NBT = NBT + 1
IBNDT (NBT) = NODE2
ITYPET(NBT) = 1
BVT (NBT) = 0.D0
C
NBT = NBT + 1
IBNDT (NBT) = NODE3
ITYPET(NBT) = 1
BVT (NBT) = 0.D0
C
IF ( I .NE. NEY ) THEN
NBT = NBT + 1
IBNDT (NBT) = NODE4
ITYPET(NBT) = 1
BVT (NBT) = 0.D0
END IF
END DO
C=======================================================================
C MAKING DATA FILES
C---------- PROJECT
C---------- PARAMETERS
OPEN ( 1, FILE=PROJECT, STATUS = 'UNKNOWN' )
WRITE (1,*) EXX, EXY, EYY
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)
END DO
C---------- DIRICHLET TYPE BOUNDARY CONDITIONS
WRITE (1,*) NBT
DO I = 1 , NBT
WRITE (1,*) IBNDT(I), ITYPET(I), BVT(I)
END DO
CLOSE (1)
C=======================================================================
C---------- ELEMENT DRAWING
OPEN ( 1, FILE='ELEMENT12.OUT', STATUS = 'UNKNOWN')
DO I = 1, NE
DO J = 1, 4
WRITE (1,*) XCOORD(NODEX(I,J)), YCOORD(NODEX(I,J))
END DO
WRITE (1,*)
DO J = 4, 7
WRITE (1,*) XCOORD(NODEX(I,J)), YCOORD(NODEX(I,J))
END DO
WRITE (1,*)
DO J = 7, 10
WRITE (1,*) XCOORD(NODEX(I,J)), YCOORD(NODEX(I,J))
END DO
WRITE (1,*)
DO J = 10, 12
WRITE (1,*) XCOORD(NODEX(I,J)), YCOORD(NODEX(I,J))
END DO
WRITE (1,*) XCOORD(NODEX(I,1)), YCOORD(NODEX(I,1))
WRITE (1,*)
END DO
CLOSE (1)
C=======================================================================
C------ CREATION OF PARAMETER FILE TO BE USED IN INCLUDE STATEMENT
CALL BANDWID ( MXE, ND, NE, NODEX, NBW )
OPEN ( 1, FILE='PARAM.DAT', STATUS='UNKNOWN' )
WRITE (1,*) ' PARAMETER ( ND=',ND,', INTEPT=',INTEPT,' )'
WRITE (1,*) ' PARAMETER ( MXE=',NE,', MXN=',NNODE,
* ', MXB=',NBT, ', MXW=',NBW, ' )'
CLOSE (1)
STOP 'NORMAL TERMINATION'
END
C
C
SUBROUTINE BANDWID ( MXE , ND , NE , NODEX , NBW )
DIMENSION NODEX(MXE,ND)
C------- RETURN VALUE: NBW
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
C
C
FUNCTION F (X,A)
IMPLICIT REAL*8 ( A-H , O-Z )
F = X**A
RETURN
END