PROGRAM SETFEM8
C=======================================================================
C DATA GENERATING PROGRAM FOR FEM4Q.FOR
C DOMAIN : TLX BY TLY
C BOUNDARY CONDITIONS
C ON FACES OF -X AND -Y: T(X,Y)=0.
C ON FACE OF +X: DT/DN = 0
C ON FACE OF +Y: T(X,Y) = SIN ( 0.25*PI*X )
C EIJI FUKUMORI NOVEMBER 28, 1995
C NUMBERING OF NODEAL POINTS
C FOR CASE OF NEX=2 & NEY=2 17-----18-----19-----20-----21
C | | |
C 14 15 16
C | | |
C 9-----10-----11-----12-----13
C | | |
C 6 7 8
C | | |
C 1------2------3------4------5
C NOEDX OF 1ST ELEMENT = 1, 2, 3, 7, 11, 10, 9, 6
C FEB. 13, 2013
C=======================================================================
IMPLICIT REAL*8 ( A-H , O-Z )
PARAMETER ( ND=8, MXE=2200, MXN=2400, MXB=1300, INTEPT=3 )
PARAMETER ( TLX = 2.D0, TLY = 2.D0 )
C=======================================================================
DIMENSION NODEX(MXE,ND),XCOORD(MXN),YCOORD(MXN),IBNDT(MXB),
* BVT(MXB), ITYPET(MXB)
CHARACTER PROJECT*14,EXFILE*3
LOGICAL YES
C=======================================================================
DATA PROJECT / 'FEM08INPUT.DAT' /
C=======================================================================
C NEY: NEMBER OF VERTICAL ELEMENTS (NUMBER OF NODES: NEY+1)
C NEX: NEMBER OF HORIZONTAL ELEMENTS (NUMBER OF NODES: NEX+1)
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=======================================================================
DX = TLX / NEX
DY = TLY / NEY
EXX = 1.D0
EXY = 0.D0
EYY = 1.D0
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)*(NEX*3+2) + 1 + (J-1)*2
NODEX(NE,2) = NODEX(NE,1) + 1
NODEX(NE,3) = NODEX(NE,2) + 1
NODEX(NE,4) = NODEX(NE,3) + (NEX*2 + 1 - J )
NODEX(NE,5) = NODEX(NE,3) + (NEX*3 + 2 )
NODEX(NE,6) = NODEX(NE,5) - 1
NODEX(NE,7) = NODEX(NE,6) - 1
NODEX(NE,8) = NODEX(NE,4) - 1
END DO
END DO
C=======================================================================
C NODAL COORDINATE CREATION
NNODE = 0
NCENTER = NEX + 1
DO I = 1 , NEY
DO J = 1 , NEX*2+1
NNODE = NNODE + 1
IF ( NNODE .GT. MXN ) STOP 'NNODE > MXN'
XCOORD(NNODE) = DX/2.D0*(J-1)
YCOORD(NNODE) = DY *(I-1)
END DO
DO J = 1 , NCENTER
NNODE = NNODE + 1
IF ( NNODE .GT. MXN ) STOP 'NNODE > MXN'
XCOORD(NNODE) = DX*(J-1)
YCOORD(NNODE) = DY*(I-1) + DY/2.D0
END DO
END DO
C-------- TOP LAYER -----------
DO J = 1 , NEX*2+1
NNODE = NNODE + 1
IF ( NNODE .GT. MXN ) STOP 'NNODE > MXN'
XCOORD(NNODE) = DX/2.D0*(J-1)
YCOORD(NNODE) = DY *NEY
END DO
C=======================================================================
PI = 4.D0* DATAN (1.D0)
C BOUNDARY CONDITIONS
C--------- MOMENTUM EQUATIONS AND HEAT EQUATION
NBT = 0
C--------- FACE OF +Y AND -Y
DO J = 1 , NEX*2+1
C--------- -Y
NBT = NBT + 1
IBNDT(NBT) = J
ITYPET(NBT) = 1
BVT(NBT) = 0.D0
C--------- +Y
NBT = NBT + 1
IBNDT(NBT) = NEY*(NEX*3 + 2) + J
ITYPET(NBT) = 1
BVT(NBT) = DSIN ( 0.25D0*PI*XCOORD(IBNDT(NBT)))
END DO
C--------- FACE OF +X AND -X
C -X
DO I = 1 , NEY
NBT = NBT + 1
IBNDT (NBT) = (I-1)*(NEX*3+2) + NEX*2 + 1 + 1
ITYPET(NBT) = 1
BVT(NBT ) = 0.D0
NBT = NBT + 1
IBNDT (NBT) = I*(NEX*3+2) + 1
ITYPET(NBT) = 1
BVT(NBT ) = 0.D0
END DO
NBT = NBT - 1
C=======================================================================
WRITE (*,*) ' NUMBER OF ELEMENTS (NE) = ',NE
WRITE (*,*) ' NUMBER OF NODAL POINTS (NNODE) = ',NNODE
C=======================================================================
C DATA FILE INQUIRY
EXFILE = 'NEW'
INQUIRE ( FILE = PROJECT, EXIST = YES )
IF ( YES ) EXFILE='OLD'
C=======================================================================
C MAKING DATA FILES
C---------- PROJECT
OPEN ( 1, FILE=PROJECT, STATUS = EXFILE )
WRITE (1,*) EXX, EXY, EYY
WRITE (1,*) NE
DO IEL = 1 , NE
WRITE (1,*) IEL,(NODEX(IEL,J),J=1,ND)
END DO
WRITE (1,*) NNODE
DO NODE = 1 , NNODE
WRITE (1,*) NODE,XCOORD(NODE),YCOORD(NODE)
END DO
WRITE (1,*) NBT
DO I = 1 , NBT
WRITE (1,*) IBNDT(I), ITYPET(I), BVT(I)
END DO
CLOSE (1)
C=======================================================================
OPEN ( 1, FILE='ELEMENT8.OUT', STATUS = 'UNKNOWN' )
DO IEL = 1 , NE
WRITE(1,*) XCOORD(NODEX(IEL,1)),YCOORD(NODEX(IEL,1))
WRITE(1,*) XCOORD(NODEX(IEL,2)),YCOORD(NODEX(IEL,2))
WRITE(1,*) XCOORD(NODEX(IEL,3)),YCOORD(NODEX(IEL,3))
WRITE(1,*)
WRITE(1,*) XCOORD(NODEX(IEL,3)),YCOORD(NODEX(IEL,3))
WRITE(1,*) XCOORD(NODEX(IEL,4)),YCOORD(NODEX(IEL,4))
WRITE(1,*) XCOORD(NODEX(IEL,5)),YCOORD(NODEX(IEL,5))
WRITE(1,*)
WRITE(1,*) XCOORD(NODEX(IEL,5)),YCOORD(NODEX(IEL,5))
WRITE(1,*) XCOORD(NODEX(IEL,6)),YCOORD(NODEX(IEL,6))
WRITE(1,*) XCOORD(NODEX(IEL,7)),YCOORD(NODEX(IEL,7))
WRITE(1,*)
WRITE(1,*) XCOORD(NODEX(IEL,7)),YCOORD(NODEX(IEL,7))
WRITE(1,*) XCOORD(NODEX(IEL,8)),YCOORD(NODEX(IEL,8))
WRITE(1,*) XCOORD(NODEX(IEL,1)),YCOORD(NODEX(IEL,1))
WRITE(1,*)
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=',NBT,')'
CLOSE (1)
C=======================================================================
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