戻る　
　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　
C
PROGRAM MAIN
C ****************************************************************** C
C C
C Fluid Analysis Software by Finite Element Method C
C C
C ( Triangular Linear Element ) C
C C
C F 2 D - F E T - FIG C
C C
C ( Version : 3.0 ) C
C C
C Copyright : Yasuhiro MATSUDA C
C C
C ****************************************************************** C
C
IMPLICIT REAL*8 (A-H, O-Z)
C
INTEGER*2 ICR, CPRT, IRWS, CON, NODE
C
C ----- ATTENTION ( 17 x 17 x 4 ) ---
C
DIMENSION X ( 113), Y ( 113), UN( 113), VN( 113), PSI( 113),
1 VOR( 113)
C
DIMENSION NODE( 196,3), IRWS( 616)
C
COMMON /CNST/ KPLT, ICR, CPRT
COMMON /TIMC/ DT, LPLM, LOP1, LOP, ALP
COMMON /CMPD/ KSY( 200), NY, KSX( 200), NX
C
C ----- Attention --- < == 500 ---
C
COMMON /VRPT/ VARV( 500), VRUV( 500), VARP( 500)
C
OPEN ( 10, FILE='F2D_FET_FIG.DAT' )
OPEN ( 8, FILE='GH100.DAT' ) ! Re = 100
C
READ(10,2040) CON, LOP, DT, RE, NNP, NEL
2040 FORMAT( I2, I4, F6.3, F8.1, 2I4)
C
C ----------------------------------------------
C
CALL PREPS ( PSI, VOR, UN, VN, X, Y, NODE, IRWS, NNP, NEL )
C
C ----------------------------------------------
C
IF ( KPLT.NE.1 ) GO TO 900
C
C ----------------------------------------------
C
CALL OTLIN ( CPRT, IRWS, X, Y, NNP )
C
C ----------------------------------------------
C
CALL PLOTN ( X, Y, UN, VN, NNP, CON )
C
C ----------------------------------------------
C
CALL PLTVR ( VRUV, LOP, 1 )
C
CALL PLTVR ( VARP, LOP, 4 )
C
CALL PLTVR ( VARV, LOP, 3 )
C
C --------------------------------------------------------------------
C
IF ( CON.EQ.0) CALL PLCNT ( X, Y, NODE, PSI, NNP, NEL, CON )
C
GO TO 990
C --------------------------------------------------------------------
C
900 CONTINUE
C
C ----------------------------------------
C
CALL CMPLT ( UN, VN, X, Y, NNP )
C
C ----------------------------------------
C
990 CONTINUE
C
CLOSE ( 8)
CLOSE ( 10)
C
STOP
END
C **********************************************************************
C
SUBROUTINE PREPS ( PSI, VOR, UN, VN, X, Y, NODE, IRWS, NNP, NEL )
C
IMPLICIT REAL*8 (A-H, O-Z)
C
INTEGER*2 ICR, CPRT, NODE, IRWS
C
DIMENSION NODE( NEL,3), IRWS( 616)
C
DIMENSION PSI( NNP), VOR( NNP), UN( NNP), VN( NNP), X( NNP),
1 Y ( NNP)
C
COMMON /CNST/ KPLT, ICR, CPRT
COMMON /TIMC/ DT, LPLM, LOP1, LOP, ALP
COMMON /VRPT/ VARV( 500), VRUV( 500), VARP( 500)
COMMON /CMPD/ KSY( 200), NY, KSX( 200), NX

C
C --------------------
C
CALL TTLE
C
C --------------------
C
ICR = 15
C
READ(10,2050) ( X( I), Y( I), I= 1, NNP )
2050 FORMAT(10F8.4)
C
READ(10,2060) ( ( NODE( I,J), I=1, NEL ), J=1, 3 )
2060 FORMAT(20I4)
READ(10,2060) CPRT, ( IRWS( I), I=1, 616 )
C
C ----- At the Last LOP ---
C
READ(10,2050) ( VRUV( I), VARP( I), I = 1, LOP )
READ(10,2050) ( VARV( I), I = 1, LOP )
C
READ(10,2050) ( UN ( I), VN ( I), I= 1, NNP )
READ(10,2050) ( PSI( I), VOR( I), I= 1, NNP )
C
C ----- To set Node numbers for a comparison with Giha's Data ---
C
NY = 0
C
DO 100 I = 1, NNP
C
IF ( DABS( X( I)-0.5D0 ).LE.0.01D0 ) THEN
NY = NY + 1 ! NY < = 200
KSY( NY) = I
END IF
C
100 CONTINUE
C
C ------------
C
NX = 0
C
DO 200 I = 1, NNP
C
IF ( DABS( Y( I)-0.5D0 ).LE.0.01D0 ) THEN
NX = NX + 1 ! NX < = 200
KSX( NX) = I
END IF
C
200 CONTINUE
C
WRITE(6,1234) NX, NY
1234 FORMAT(/,' *** NX =',I3,' NY =',I3)
C
C ----------------------------
C
WRITE(6,*) ' '
WRITE(6,*) ' * KPLT = 1 ( Meshes / Str.Fn./Vel.Dis. ) or 2 ( Compa
1rison with Ref. )'
READ(5,*) KPLT
WRITE(6,*) ' '
WRITE(6,*) ' * KPLT =', KPLT
C
CALL VGA@
C
RETURN
END
C *********************************************************************
C
SUBROUTINE CMPLT ( UN, VN, X, Y, NNP )
C
IMPLICIT REAL*8(A-H,O-Z)
C
INTEGER*2 ICR, IX1, IY1, IX2, IY2
C
DIMENSION UN( NNP), VN( NNP), X( NNP), Y( NNP)
C
COMMON /TIMC/ DT, LPLM, LOP1, LOP, ALP
C
NPL = 0
C ---------------
C
IX1 = 0
IY1 = 0
IX2 = 639
IY2 = 479
ICR = 7
C
CALL CLEAR_SCREEN_AREA@ ( IX1, IY1, IX2, IY2, ICR )
C
KEY = 0
C
CALL REFCMP ( X, Y, UN, VN, NNP, KEY, NPL )
C
KEY = KEY + 1
C
RETURN
END
C **********************************************************************
C
SUBROUTINE REFCMP ( X, Y, U, V, NNP, KEY, NPL )
C
IMPLICIT REAL*8(A-H,O-Z)
C
DIMENSION X( NNP), Y( NNP), U( NNP), V( NNP)
C
DIMENSION GU( 17), GV( 17), GUY( 17), GVX( 17)
C
COMMON /CMPD/ KSY( 200), NY, KSX( 200), NX
C
INTEGER*2 IX1, IY1, IX2, IY2, ICR, IICR
C
CHARACTER*80 STR0
C
DATA STR0 /'** Results of F2D-FET * ( Velocity Distribution compa
1red with Ghia''s Data )'/
C
IF ( KEY.EQ.0) ICR = 11
IF ( KEY.EQ.1) ICR = 13
IF ( KEY.EQ.2) ICR = 14
C
IICR = 16
C
CALL SET_TEXT_ATTRIBUTE@ ( 3, 1., 0., 0. )
C
CALL DRAW_TEXT@ ( STR0, 10, 5, IICR )
C
C ----- Ghia's Data --- ( Re = 100 or 400 or 1000 )---
C
DO 100 I = 1, 17
C
READ(8,1000) GU( I), GUY( I), GV( I), GVX( I)
1000 FORMAT(4F10.5)
C
100 CONTINUE
C
REWIND ( 8)
C
C ----- Attention : For Velocity direction change ---
C
DO 200 I = 1, 17
C
GU ( I) = - GU( I)
GVX( I) = 1.D0 - GVX( I)
C
200 CONTINUE
C ----------------------------
C
X1 = 345
Y1 = 50
X0 = 20
Y0 = 50
C
SAL = 260.D0
X2 = SAL* 0.5D0
Y2 = SAL* 0.5D0
C
C ----- U by Ghia ---
C
IX1 = X0
IY1 = Y0
IX2 = X0 + SAL
IY2 = Y0 + SAL
C
CALL RECTANGLE@ ( IX1, IY1, IX2, IY2, IICR )
C
IX1 = X0 + X2
IY1 = Y0 + SAL
IX2 = X0 + X2
IY2 = Y0
C
CALL RECTANGLE@ ( IX1, IY1, IX2, IY2, IICR )
C
C ----------------------------------------------------
DO 300 I = 1, 16
C
IC = I + 1
C
IX1 = X0 + X2 + GU ( I)* X2
IY1 = Y0 + SAL - GUY( I)* SAL
IX2 = X0 + X2 + GU ( IC)* X2
IY2 = Y0 + SAL - GUY( IC)* SAL
C
CALL DRAW_LINE@ ( IX1, IY1, IX2, IY2, IICR )
C
300 CONTINUE
C ----------------------------------------------------
C
IX1 = 40
IY1 = 80 + NPL*20
C
IX2 = 80
IY2 = IY1
C
CALL DRAW_LINE@ ( IX1, IY1, IX2, IY2, IICR )
CALL DRAW_TEXT@ ( ' Ref.', IX2+5, IY1-7, IICR )
C
NPL = NPL + 1
C -------------------
C
C ***** Computed U *******************************
C
DO 400 J = 1, NY
C
NN = KSY( J)
C
IX1 = X0 + X2 + U( NN)* X2
IY1 = Y0 + SAL - Y( NN)* SAL
C
IX2 = 2
IY2 = 2
C
CALL ELLIPSE@ ( IX1, IY1, IX2, IY2, ICR )
C
400 CONTINUE
C ************************************************
C
C ----- V by Ghia ---
C
IX1 = 345
IY1 = 50
IX2 = X1 + SAL
IY2 = Y1 + SAL
C
CALL RECTANGLE@ ( IX1, IY1, IX2, IY2, IICR )
C
IX1 = X1
IY1 = Y1 + Y2
IX2 = X1 + SAL
IY2 = Y1 + Y2
C
CALL RECTANGLE@ ( IX1, IY1, IX2, IY2, IICR )
C
DO 500 I = 1, 16
C
IC = I + 1
C
IX1 = X1 + ( 1 - GVX( I))* SAL
IY1 = Y1 + Y2 + GV ( I) * Y2
C
IX2 = X1 + ( 1 - GVX( IC))* SAL
IY2 = Y1 + Y2 + GV ( IC) * Y2
C
CALL DRAW_LINE@ ( IX1, IY1, IX2, IY2, IICR )
C
500 CONTINUE
C
C ***** Computed V ***********************************
C
C ----------------------------
DO 600 J = 1, NX
C
NN = KSX( J)
C
IX1 = X1 + X( NN)* SAL
IY1 = Y1 + Y2 - V( NN)* Y2
C
IX2 = 2
IY2 = 2
C
CALL ELLIPSE@ ( IX1, IY1, IX2, IY2, ICR )
C
600 CONTINUE
C ----------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE TTLE
C
IMPLICIT REAL*8 ( A-H, O-Z)
C
WRITE(6,*) ' '
WRITE(6,*) '****************************************************'
WRITE(6,*) ' '
WRITE(6,*) ' Fluid Analysis Software by FEM '
WRITE(6,*) ' '
WRITE(6,*) ' F I G - F 2 D - F E T '
WRITE(6,*) ' '
WRITE(6,*) ' ( Modified Galerkin Method - Triangular Linear ) '
WRITE(6,*) ' '
WRITE(6,*) ' ( Version: 3.0 ) '
WRITE(6,*) ' '
WRITE(6,*) ' Copyright 2011 : Yasuhiro MATSUDA '
WRITE(6,*) ' '
WRITE(6,*) '****************************************************'
C
RETURN
END
C ***********************************************************************
C
SUBROUTINE PLCNT ( X, Y, NODE, TEP, NNP, NEL, CON )
C
IMPLICIT REAL*8 (A-H, O-Z)
C
C ----- To plot contour lines for Triangular Meshes ---
C
C ( Linear interpolation ) ---
C
INTEGER*2 NODE, LC, CON
C
INTEGER*2 ICR, IWD, IX0, IY0, IX1, IY1, IX2, IY2, IX3, IY3
C
DIMENSION X( NNP), Y( NNP), NODE(NEL,3), TEP( NNP), PX( 3),
1 PY( 3), LC( 9), HVN( 20)
C
ICR = 15
C
IWD = 160
C
IX0 = 280
IY0 = 70
IX1 = IX0 + IWD
IY1 = IY0 + IWD
C
CALL RECTANGLE@ ( IX0, IY0, IX1, IY1, ICR )
C
LC( 1) = 1
LC( 2) = 5
LC( 3) = 2
LC( 4) = 6
LC( 5) = 4
C
LC( 6) = 7
LC( 7) = 7
LC( 8) = 7
LC( 9) = 7
C
NM = 10
C
C ---------------------------------------------------------------
C
TMIN = 100.D0
TMAX = -100.D0
C
DO 30 I = 1, NNP
C
IF ( TEP( I).GE.TMAX ) TMAX = TEP( I)
IF ( TEP( I).LE.TMIN ) TMIN = TEP( I)
C
30 CONTINUE
C
WRITE(6,*) ' '
WRITE(6,1356) TMAX, TMIN
1356 FORMAT(' * Cmax.=',F7.3,' * Cmin.=',F7.3)
C
C ----------------------------
C
RTD = ( TMAX - TMIN ) / DFLOAT( NM)
DO 50 I = 1, NM+1
HVN( I) = TMIN + RTD* DFLOAT( I-1)
IF ( HVN( I).LT.0.D0 ) NMZ = I
50 CONTINUE
C
*** WRITE(6,*) ' * HVN=', ( HVN( I), I=1, NM+1 )
C ----------------------------
C
ICR = 11
C
C ----- Plot contour lines ---
C
PXV = 0.D0
PYV = 0.D0
C
C ------------------------------------------
DO 100 L = 1, NEL
C
C ---------------------------------
DO 200 IH = 1, NM
C
HV = HVN( IH)
C
ICL = 0
C
C ----------------------------
DO 300 J = 1, 3
C
I1 = NODE( L, J)
C
IF ( J.NE.3) THEN
C
I2 = NODE( L,J+1)
C
ELSE
I2 = NODE( L, 1)
C
END IF
C
IF ( DABS( TEP( I1) - TEP( I2) ).LE. 0.00001D0 ) THEN
IF ( DABS( TEP( I1) - HV ) .LE. 0.00001D0 ) THEN
C
ICL = ICL + 1
C
PX( ICL) = X( I1)* DFLOAT( IWD)
PY( ICL) = Y( I1)* DFLOAT( IWD)
END IF
C
ELSE IF (((TEP( I1).LE. HV) .AND. (HV .LT. TEP( I2)))
1 .OR. ((TEP( I1).GE. HV) .AND. (HV .GT. TEP( I2)))) THEN
C
ICL = ICL + 1
RT = (HV - TEP( I1)) / (TEP( I2) - TEP( I1))
C
PX( ICL) = ( X( I1) + ( X( I2) - X( I1) )* RT )* DFLOAT( IWD)
PY( ICL) = ( Y( I1) + ( Y( I2) - Y( I1) )* RT )* DFLOAT( IWD)
C
END IF
C
300 CONTINUE
C ----------------------------
C
IF ( ICL .GE. 2 ) THEN
C
C ----------------------------
C
DO 400 J = 1, ICL
C
JM = J
D2 = ( PXV - PX( J) )** 2 + ( PYV - PY( J) )** 2
C
C -----------------------
C
DO 500 JJ = J+1, ICL
C
D2V = ( PXV - PX( JJ) )** 2 + ( PYV - PY( JJ) )** 2

IF ( D2V .LT. D2 ) THEN
C
JM = JJ
D2 = D2V
C
END IF
C
500 CONTINUE
C -----------------------
C
PXV = PX( JM)
PYV = PY( JM)
C
PX( JM) = PX( J)
PY( JM) = PY( J)
C
PX( J) = PXV
PY( J) = PYV
C
400 CONTINUE
C ----------------------------
C
PXV = PX( ICL)
PYV = PY( ICL)
C
END IF
C ----------------------------
IF ( ICL.GT.0) THEN
C
DO 600 J = 2, ICL
C
IX1 = IX0 + IDNINT( PX( J-1))
IY1 = IY0 + IWD - IDNINT( PY( J-1))
C
IX2 = IX0 + IDNINT( PX( J))
IY2 = IY0 + IWD - IDNINT( PY( J))
C
CALL DRAW_LINE@ ( IX1, IY1, IX2, IY2, ICR )
C
600 CONTINUE
C
IF ( ICL.EQ.3) THEN
IX3 = IX0 + IDNINT( PX( 1))
IY3 = IY0 + IWD - IDNINT( PY( 1))
C
CALL DRAW_LINE@ ( IX2, IY2, IX3, IY3, ICR )
C
END IF
C
END IF
C ----------------------------
C
200 CONTINUE
C ---------------------------------
C
100 CONTINUE
C ------------------------------------------
C
RETURN
END
C ***********************************************************************
C
SUBROUTINE PLOTN ( X, Y, UN, VN, NNP, CON )
C
IMPLICIT REAL*8 (A-H, O-Z)
C
INTEGER*2 CON, ICR, IWD, IX0, IY0, IX1, IY1, IX2, IY2, GX, GY
C
DIMENSION X( NNP), Y( NNP), UN( NNP), VN( NNP)
C
COMMON /TIMC/ DT, LPLM, LOP1, LOP, ALP
C
SCLX = 240.D0
SCLY = 240.D0
SCL = 50.D0
C
ICR = 15
C
IX0 = 460
IY0 = 70
C
IWD = 160
c
TL = DT* DFLOAT( LOP)
C
IX1 = IX0 + IWD
IY1 = IY0 + IWD
C
CALL RECTANGLE@ ( IX0, IY0, IX1, IY1, ICR )
C
ICR = 11
C
DO 100 I = 1, NNP
C
GX = IX0 + X( I)* DFLOAT( IWD)
GY = IY0 + IWD - Y( I)* DFLOAT( IWD)
C
IF ( GX.EQ.IX0 ) GO TO 100
IF ( GY.EQ.IY0 ) GO TO 100
IF ( GX.EQ.IX0+IWD ) GO TO 100
IF ( GY.EQ.IY0+IWD ) GO TO 100
C
IX1 = GX
IY1 = GY
IX2 = GX + IDNINT( UN( I)* SCL )
IY2 = GY - IDNINT( VN( I)* SCL )
C
CALL DRAW_LINE@ ( IX1, IY1, IX2, IY2, ICR )
C
100 CONTINUE
C
RETURN
END
C **********************************************************************
C
SUBROUTINE OTLIN ( NLIN, ILINE, X, Y, NNP )
C
IMPLICIT REAL*8 (A-H, O-Z)
C
C ----- MESH PLOT ---
C
INTEGER*2 ILINE, NLIN, ICR, IWD, IX0, IY0, IX1, IY1, IX2, IY2
C
DIMENSION ILINE( NLIN), X( NNP), Y( NNP)
C
ICR = 15
C
IX0 = 100
IY0 = 70
IWD = 160
C
DO 100 I = 1, NLIN
C
N = ILINE( ( I-1)* 2 + 1 )

IX1 = IDNINT( X( N)* DFLOAT( IWD) ) + IX0
IY1 = IDNINT( Y( N)* DFLOAT( IWD) ) + IY0
C
N = ILINE( I*2)
C
IX2 = IDNINT( X( N)* DFLOAT( IWD) ) + IX0
IY2 = IDNINT( Y( N)* DFLOAT( IWD) ) + IY0
C
CALL DRAW_LINE@ ( IX1, IY1, IX2, IY2, ICR )
C
100 CONTINUE
C
RETURN
END
C **********************************************************************
C
SUBROUTINE PLTVR ( RAX, NPL, KKK )
C
IMPLICIT REAL*8(A-H,O-Z)
C
INTEGER*2 ICR, IX0, IY0
C
C ----- Attention NPL .LE. 500 ---
C
DIMENSION RAX( NPL)
C
INTEGER*2 IVAR( 500), IXX( 500)
C
CHARACTER*80 STR3
CHARACTER*6 STR1( 3)
CHARACTER*3 STR2( 4)
C
DATA STR3 /'* F2D-FET( SV ) ( Max.Var_of Str. funct., u-Velocity
1 & Vorticity(%) )'/
DATA STR1/'0 ','NPL/2',' NPL'/
DATA STR2/'200','150','100',' 50'/
C
C ----- Attention --- Max. Value ---
C
IF ( KKK.GE.2) GO TO 1234
C
ICR = 15
C
CALL SET_TEXT_ATTRIBUTE@ ( 3, 1., 0., 0. )
C
CALL DRAW_TEXT@ ( STR3, 70, 15, ICR )
C
ICR = 15
C
CALL DRAW_TEXT@ ( STR1( 1), 47, 450, ICR )
CALL DRAW_TEXT@ ( STR1( 2), 305, 450, ICR )
CALL DRAW_TEXT@ ( STR1( 3), 585, 450, ICR )
C
IF ( KKK.NE.1 ) GO TO 133
C
CALL DRAW_TEXT@ ( '1000', 15, 35, ICR )
CALL DRAW_TEXT@ ( ' 100', 15, 115, ICR )
CALL DRAW_TEXT@ ( ' 10', 15, 195, ICR )
CALL DRAW_TEXT@ ( ' 1', 15, 275, ICR )
CALL DRAW_TEXT@ ( ' 0.1', 15, 355, ICR )
CALL DRAW_TEXT@ ( '0.01', 15, 435, ICR )
C
CALL DRAW_TEXT@ ( '(%)', 15, 240, ICR )
C
133 CONTINUE
C
IX0 = 50
IY0 = 40
C
CALL RECTANGLE@ ( IX0, IY0, IX0+580, IY0+400, 15 )
C
1234 CONTINUE
C
IX0 = 50
IY0 = 40
C -------------
C
RIWD = 580.D0 / DFLOAT( NPL)
C =================================
C
IF ( KKK.EQ.1) ICR = 15
IF ( KKK.EQ.2) ICR = 12
IF ( KKK.EQ.3) ICR = 9
IF ( KKK.EQ.4) ICR = 11
C
C --------------------------------------
C
DO 100 I = 1, NPL
IF ( RAX( I).LE.0.D0 ) GO TO 100
RAX( I) = DLOG10( RAX( I))
100 CONTINUE
C
C --------------------------------------
C
DO 200 I = 1, NPL
C
IXX( I) = IX0 + IDNINT( DFLOAT( I)* RIWD )
C
IVAR( I) = 280 - IDNINT( 80.D0* RAX( I) )
C
200 CONTINUE
C --------------------------------------
C
CALL POLYLINE@ ( IXX, IVAR, NPL, ICR )
C
*** WRITE(6,*) ' * RAX=', RAX(NPL),' IVAR(NPL)=', IVAR( NPL)
C
RETURN
END
C **********************************************************************