PROGRAM SETSTAX12 C======================================================================= C DATA GENERATING PROGRAM FOR STATICX.FOR C PROJECT NAME: C BOUNDARY CONDITIONS C LEFT END: FIXED (EMULATION OF SYMENTRIC BEAM CEMTER) C RIGHT END: MOMENT IS APPLIED C ALL OTHERS: NO NORMAL FOECE, NO SHEAR FORCE C DOMAIN: RECTANGULAR C ELEMENT: 12-NODED ISOPARAMETRIC ELEMENT C DOMAIN DISCRETIZATION: ONE-SIZE ELEMENTS C EIJI FUKUMORI APRIL 13, 2010 C======================================================================= IMPLICIT REAL*8 ( A-H , O-Z ) PARAMETER ( ND=12, INTEPT=4, MXE=22000, MXN=23000, MXB=2000 ) PARAMETER (YOUNG=69.D6,POISSON=0.3D0, HEIGHT=2.D0,WIDTH=20.D0, * TLX = WIDTH, TLY = HEIGHT, NEY = 16, NEX = 100, * DX = TLX / NEX, DY = TLY / NEY) C======================================================================= DIMENSION NODEX(MXE,ND), XCOORD(MXN), YCOORD(MXN), * IBNDFX(MXB),IBNDFY(MXB),BVX(MXB),BVY(MXB), * IBFORCEX(MXB), IBFORCEY(MXB), BVFORCEX(MXB), BVFORCEY(MXB) DIMENSION NEUTRAL(MXN) CHARACTER INPFILE*12,EXFILE*3 LOGICAL YES C======================================================================= DATA INPFILE / 'STATIC12.DAT' / C======================================================================= F(X) = X G(X) = X C======================================================================= C NEY: NEMBER OF VERTICAL ELEMENTS (NUMBER OF NODES: NEY+1) C NEX: NEMBER OF HORIZONTAL ELEMENTS (NUMBER OF NODES: NEX+1) C HEIGHT: HEIGHT OF DRIVEN CAVITY FLOW DOMAIN C======================================================================= WRITE (*,*)' YOUNG MODULUS = ',YOUNG WRITE (*,*)' POISSON RATIO = ',POISSON C======================================================================= C ELEMENT CREATION NDYEG = NEY*3+1 NDYCT = NEY + 1 NE = 0 DO J = 1 , NEX DO I = 1 , NEY NE = NE + 1 IF ( NE .GT. MXE ) STOP 'NE > MXE' NODEX(NE, 1) = (NDYEG+NDYCT*2)*(J-1) + (I-1)*3 + 1 NODEX(NE, 2) = NODEX(NE, 1) + (NEY-I+1)*3 + I NODEX(NE, 3) = NODEX(NE, 2) + NEY + 1 NODEX(NE, 7) = NODEX(NE, 1) + NDYEG + 2*NDYCT + 3 NODEX(NE, 6) = NODEX(NE, 7) - 1 NODEX(NE, 5) = NODEX(NE, 6) - 1 NODEX(NE, 4) = NODEX(NE, 5) - 1 NODEX(NE, 8) = NODEX(NE, 3) + 1 NODEX(NE, 9) = NODEX(NE, 2) + 1 NODEX(NE,12) = NODEX(NE, 1) + 1 NODEX(NE,11) = NODEX(NE,12) + 1 NODEX(NE,10) = NODEX(NE,11) + 1 END DO END DO C======================================================================= C NODAL COORDINATE CREATION DYEG = DY/3.D0 DXEG = DX/3.D0 NNODE = 0 C DO I = 1 , NEX C DO J = 1 , NDYEG NNODE = NNODE + 1 IF ( NNODE .GT. MXN ) STOP 'NNODE > MXN' YCOORD(NNODE) = DYEG * (J-1) - TLY/2.D0 XCOORD(NNODE) = (I-1)* DX END DO C DO J = 1 , NDYCT NNODE = NNODE + 1 IF ( NNODE .GT. MXN ) STOP 'NNODE > MXN' YCOORD(NNODE) = DY * (J-1) - TLY/2.D0 XCOORD(NNODE) = (I-1)* DX + DXEG END DO C DO J = 1 , NDYCT NNODE = NNODE + 1 IF ( NNODE .GT. MXN ) STOP 'NNODE > MXN' YCOORD(NNODE) = DY * (J-1) - TLY/2.D0 XCOORD(NNODE) = (I-1)* DX + DXEG + DXEG END DO END DO C----------- LAST EADGE -------- DO J = 1 , NDYEG NNODE = NNODE + 1 IF ( NNODE .GT. MXN ) STOP 'NNODE > MXN' YCOORD(NNODE) = DYEG * (J-1) - TLY/2.D0 XCOORD(NNODE) = NEX * DX END DO C======================================================================= C BOUNDARY CONDITIONS C==== FIRST KIND C--------- NAVIER EQUATIONS NBFX = 0 NBFY = 0 NFORCEX = 0 NFORCEY = 0 C--------- X DIRECTION ON FACE -X U0 = 0.D0 DO J = 1 , NDYEG NBFX = NBFX + 1 IBNDFX(NBFX) = J BVX(NBFX) = U0 END DO C======================================================== C--------- Y DIRECTION ON FACE -X NODE AT CENTER NODE C======================================================== V0 = 0.D0 NBFY = NBFY + 1 IBNDFY(NBFY) = NEY/2*3 + 1 BVY(NBFY) = V0 C======================================================== C======================================================== C==== SECOND KIND C--------- FACE OF +X FORCEINX = 1000.D0 C LOWER POINT NFORCEX = NFORCEX + 1 IBFORCEX(NFORCEX) = NNODE - NDYEG + 1 BVFORCEX(NFORCEX) = FORCEINX C UPPER POINT FORCEINX = -1000.D0 NFORCEX = NFORCEX + 1 IBFORCEX(NFORCEX) = NNODE BVFORCEX(NFORCEX) = FORCEINX C======================================================================= C NODEL NUMBERS ALONG NEUTRAL LINE OF BEAM N1 = NEY/2*3 N2 = NEY/2 NNEUT = 1 NEUTRAL(NNEUT) = N1 + 1 DO I = 1 , NEX NNEUT = NNEUT + 1 NEUTRAL(NNEUT) = NEUTRAL(NNEUT-1) + N1 + N2 + 1 NNEUT = NNEUT + 1 NEUTRAL(NNEUT) = NEUTRAL(NNEUT-1) + 2*N2 + 1 NNEUT = NNEUT + 1 NEUTRAL(NNEUT) = NEUTRAL(NNEUT-1) + N2 + N1 + 1 END DO 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) 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 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 C---------- FINAL CLOSE (1) C======================================================================= C---------- ELEMENT DRAWING OPEN ( 1, FILE='ELEMENT12.DAT', 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---------- NEUTRAL NODAL NUMBERS OPEN ( 1, FILE='NEUTRAL12.DAT', STATUS = 'UNKNOWN') WRITE (1,*) NNEUT DO I = 1 , NNEUT WRITE (1,*) NEUTRAL(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,')' NB = MAX0 (NBFX, NBFY, NFORCEX, NFORCEY ) WRITE (1,*) ' PARAMETER ( MXB=',NB,' )' 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