戻る

C
PROGRAM MAIN
C
C ******************************************************************** C
C C
C Graphics for F 2 D - F E M C
C C
C ( Version 5.1 ) C
C C
C Copyright : Yasuhiro MATSUDA C
C C
C ( 2D Fluid Flow Analysis by FEM ) C
C C
C ( uvp / Quadrilateral Elements ) C
C C
C ******************************************************************** C
C
IMPLICIT REAL*8 (A-H,O-Z)
C
C ----- Attention --- Max. for 50 x 50 --- NNP: 2601, NEL: 2500 ---
C ==================
C
DIMENSION X( 2601), Y ( 2601), UN( 2601), VN( 2601),
1 P( 2601), NDE( 2500,4)
C
C --------------------------------------------------------------------
C
COMMON /CNTL/ CON
COMMON /TIMC/ DT, LPLM, LOP, ALP, RE
C
C ----- < = 3000 ---
C
COMMON /VRPL/ VRV( 3000)
C
C -----------------
C
** OPEN ( 11, FILE ='UM30_R100_F2D_FEM_PLT.OUT') ! DT = 0.04
** OPEN ( 8, FILE='GH100.DAT' ) ! Re = 100
C
C -----------------
C
** OPEN ( 11, FILE='UM30_R400_F2D_FEM_PLT.OUT') ! Re= 400 DT= 0.03
** OPEN ( 8, FILE='GH400.DAT' )
C
C -----------------
C
OPEN ( 11, FILE ='UM30_R1000_002_F2D_FEM_PLT.OUT')
C ! Re = 1000 DT = 0.02
* OPEN ( 11, FILE ='UM30_R1000_003_F2D_FEM_PLT.OUT') ! DT = 0.03
C
OPEN ( 8, FILE='GH1000.DAT' ) ! Re = 1000
C
C --------------------------------------------------------------------
C
CALL TTLE
C
C -----------------
READ(11,1000) LOP, DT, RE, NNP, NEL
1000 FORMAT( I4, F6.3, F8.1, 2I4)
C
CALL SBIP ( X, Y, UN, VN, P, NDE, NNP, NEL, KPT )
C
CALL VGA@
C
IF ( KPT.NE.1) GO TO 333
C
CALL MESF ( X, Y, NDE, Xmin, NNP, NEL )
C
CALL PLTN ( X, Y, UN, VN, NNP )
C
CALL PLVR ( VRV, LOP, 1 )
C
CALL CNTF ( X, Y, NDE, P, NNP, NEL )
C
GO TO 999
333 CONTINUE
NPL = 0
C
CALL CMPL ( UN, VN, X, Y, NNP, NPL )
C
999 CONTINUE
C
CLOSE (11)
CLOSE ( 8)
C
STOP
END
C **********************************************************************
C
SUBROUTINE CMPL ( UN, VN, X, Y, NNP, NPL )
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, LOP, ALP, RE
C
NPL = 0
C
IX1 = 0
IY1 = 0
IX2 = 639
IY2 = 479
C
ICR = 7
C
CALL CLEAR_SCREEN_AREA@( IX1, IY1, IX2, IY2, ICR )
C
KEY = 0
C
C ----------------------------------------------------------
C
CALL REFCMP ( X, Y, UN, VN, NNP, KEY, NPL )
C
C ----------------------------------------------------------
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)
DIMENSION GU( 17), GV( 17), GUY( 17), GVX( 17)
C
COMMON /TIMC/ DT, LPLM, LOP, ALP, RE
COMMON /CMPD/ KSY( 200), NY, KSX( 200), NX
C
INTEGER*2 IX1, IY1, IX2, IY2, ICR, ICR2
C
CHARACTER*80 STR0
C
DATA STR0 /'** Results of F2D-FEM * ( 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
ICR2 = 16
C
CALL SET_TEXT_ATTRIBUTE@( 3, 1., 0., 0. )
C
CALL DRAW_TEXT@( STR0, 10, 5, ICR2 )
C
C ----- Ghia's Data --- ( Re = 100 or 400 or 1000 )---
C
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 -----------------------------------------------------
C
REWIND ( 8)
C
C ----- Attention : For Velocity direction change ---
C
C -----------------------------------------------------
IF ( DABS( RE - 1000.D0 ).LE.0.001D0) THEN
C
DO 200 I = 1, 17
C
GV ( I) = - GV( I)
GVX( I) = 1.D0 - GVX( I)
C
200 CONTINUE
C
END IF
C -----------------------------------------------------
C
C -----------------------------------------------------
IF ( DABS( RE - 400.D0).LE.0.001D0) THEN
C
DO 300 I = 1, 17
C
GV ( I) = - GV( I) ! for Re= 1000
GVX( I) = 1.D0 - GVX( I)
C
300 CONTINUE
C
END IF
C -----------------------------------------------------
IF ( DABS( RE-100.D0).LE.0.001D0) THEN
C
DO 400 I = 1, 17
C
GU ( I) = - GU( I)
GVX( I) = 1.D0 - GVX( I)
C
400 CONTINUE
C
END IF
C -----------------------------------------------------
C
X1 = 345
Y1 = 50
X0 = 20
Y0 = 50
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
CALL RECTANGLE@( IX1, IY1, IX2, IY2, ICR2 )
C
IX1 = X0 + X2
IY1 = Y0 + SAL
IX2 = X0 + X2
IY2 = Y0
CALL RECTANGLE@( IX1, IY1, IX2, IY2, ICR2 )
C
C -----------------------------------------------------
DO 500 I = 1, 16
C
IC = I + 1
IX1 = X0 + X2 + GU ( I)* X2
IY1 = Y0 + SAL - GUY( I)* SAL
IX2 = X0 + X2 + GU (IC)* X2
IY2 = Y0 + SAL - GUY(IC)* SAL
CALL DRAW_LINE@( IX1, IY1, IX2, IY2, ICR2 )
C
500 CONTINUE
C -----------------------------------------------------
C
IX1 = 40
IY1 = 80 + NPL* 20
IX2 = 80
IY2 = IY1
CALL DRAW_LINE@( IX1, IY1, IX2, IY2, ICR2 )
CALL DRAW_TEXT@( ' Ref.', IX2+5, IY1-7, ICR2 )
NPL = NPL + 1
C
C ----- Computed U -------------------------------
DO 600 J = 1, NY
C
NN = KSY( J)
IX1 = X0 + X2 + U( NN)* X2
IY1 = Y0 + SAL - Y( NN)* SAL
IX2 = 2
IY2 = 2
C
CALL ELLIPSE@( IX1, IY1, IX2, IY2, ICR )
C
600 CONTINUE
C ------------------------------------------------
C
C ----- V by Ghia ---
C
IX1 = 345
IY1 = 50
IX2 = X1 + SAL
IY2 = Y1 + SAL
CALL RECTANGLE@( IX1, IY1, IX2, IY2, ICR2 )
C
IX1 = X1
IY1 = Y1 + Y2
IX2 = X1 + SAL
IY2 = Y1 + Y2
CALL RECTANGLE@( IX1, IY1, IX2, IY2, ICR2 )
C
C ------------------------------------------------
DO 700 I = 1, 16
C
IC = I + 1
IX1 = X1 + ( 1.D0 - GVX( I))* SAL
IY1 = Y1 + Y2 + GV ( I) * Y2
IX2 = X1 + ( 1.D0 - GVX( IC))* SAL
IY2 = Y1 + Y2 + GV ( IC) * Y2
CALL DRAW_LINE@( IX1, IY1, IX2, IY2, ICR2 )
C
700 CONTINUE
C ------------------------------------------------
C
C ----- Computed V -------------------------------
C
C ------------------------------------------------
DO 800 J = 1, NX
C
NN = KSX( J)
IX1 = X1 + X( NN)* SAL
IY1 = Y1 + Y2 - V( NN)* Y2
IX2 = 2
IY2 = 2
CALL ELLIPSE@( IX1, IY1, IX2, IY2, ICR )
C
800 CONTINUE
C ------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE SBIP ( X, Y, UN, VN, P, NDE, NNP, NEL, KPT )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
DIMENSION X ( NNP), Y( NNP), UN( NNP), VN( NNP), P( NNP),
1 NDE( NEL,4)
C
COMMON /CMPD/ KSY( 200), NY, KSX( 200), NX
C
C ----- Attention : Max.No_Meshes < = 200 --------
C
COMMON /VRPL/ VRV( 3000)
COMMON /TIMC/ DT, LPLM, LOP, ALP, RE
C
WRITE(6,*) ' '
WRITE(6,2000) LOP, DT, RE, NNP, NEL
2000 FORMAT( ' * LOP =',I4,' DT =',F6.3,' Re =',F8.1,' NNP =',
1 I4,' NEL =', I4 )
WRITE(6,*) ' '
WRITE(6,*) ' ** KPT = 1 ( Meshes / Vel. & Pr. Distr.) KPT= 2
1( Accuracy Check ) '
READ(5,*) KPT
WRITE(6,*) ' * KPT =', KPT
C
READ(11,1000) ( X( I), Y( I), I = 1, NNP )
1000 FORMAT(10F9.4)
READ(11,1100) ( ( NDE( I,J), J = 1, 4), I = 1, NEL )
1100 FORMAT(20I4)
READ(11,1000) ( VRV( I), I = 1, LOP )
READ(11,1000) ( UN( I), VN( I), I = 1, NNP )
READ(11,1000) ( P ( I), I = 1, NNP )
C
C ----- To set Node numbers for a comparison with Giha's Data ---
C
C ------------------------------------------------
NY = 0
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
** WRITE(6,*) ' * NY =', NY,' KSY =', ( KSY( I), I= 1, NY )
C
C ------------------------------------------------
NX = 0
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,*) ' * NX =', NX,' KSX =', ( KSX( I), I= 1, NX )
** WRITE(6,*) ' * NY =', NY,' NX =', NX
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 F.E.M. '
WRITE(6,*) ' '
WRITE(6,*) ' Graphics for F 2 D - F E M '
WRITE(6,*) ' '
WRITE(6,*) ' ( V.5.1 ) '
WRITE(6,*) ' '
WRITE(6,*) ' ( Modified Galerkin Method ) '
WRITE(6,*) ' '
WRITE(6,*) ' ( uvp / Quadrilateral Elements ) '
WRITE(6,*) ' '
WRITE(6,*) ' Copyright 2011 : Yasuhiro MATSUDA '
WRITE(6,*) ' '
WRITE(6,*) ' ****************************************************'
C
RETURN
END
C **********************************************************************
C
SUBROUTINE MESF ( X, Y, NDE, XMIN, NNP, NEL )
C
IMPLICIT REAL*8 (A-H,O-Z)

INTEGER*2 IX1, IX2, IY1, IY2, ICR
C
C ----- Important ---
C
DIMENSION X( NNP), Y( NNP), NDE( NEL,4)
C
CX = 100.D0
CY = 230.D0
C
C ------------------------------------------------
DO 100 I = 1, NEL
C
C --------------------------------------
DO 50 J = 1, 4
C
JJ = J + 1
IF ( J.EQ.4 ) JJ = 1
IX1 = CX + X( NDE( I,J ))* 160.D0
IY1 = CY - Y( NDE( I,J ))* 160.D0
IX2 = CX + X( NDE( I,JJ))* 160.D0
IY2 = CY - Y( NDE( I,JJ))* 160.D0
ICR = 14
CALL DRAW_LINE@( IX1, IY1, IX2, IY2, ICR )
C
50 CONTINUE
C --------------------------------------
C
100 CONTINUE
C ------------------------------------------------
C
Xmin = 1.D0
C
C ------------------------------------------------
DO 200 I = 2, NNP
C
XX = X( I)
AXX = DABS( X( I))
IF ( AXX.GE.0.0001D0 ) XMIN = DMIN1( XMIN,XX)
C
200 CONTINUE
C ------------------------------------------------
WRITE(6,2000) Xmin
2000 FORMAT(/,' * Xmin.=',F7.4)
C
RETURN
END
C **********************************************************************
C
SUBROUTINE PLTN ( X, Y, UN, VN, NNP )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
INTEGER*2 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, LOP, ALP, RE
C
SCL = 20.D0
IX0 = 280
IY0 = 70
IWD = 160
C
ICR = 15
TL = DT* DFLOAT( LOP)
IX1 = IX0 + IWD
IY1 = IY0 + IWD
CALL RECTANGLE@( IX0, IY0, IX1, IY1, ICR )
C
ICR = 11
C
C -----------------------------------------------------
DO 100 I = 1, NNP
C
GX = IX0 + NINT( X( I)* DFLOAT( IWD))
GY = IY0 + IWD - NINT( Y( I)* DFLOAT( IWD))
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 + NINT( UN( I)* SCL )
IY2 = GY - NINT( VN( I)* SCL )
CALL DRAW_LINE@( IX1, IY1, IX2, IY2, ICR )
C
100 CONTINUE
C -----------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE PLVR ( VRV, LOP, KEY )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
INTEGER*2 ICR, IX0, IY0
C
C ----- Attention N.LE.3000 ---
C
DIMENSION VRV( 3000)
C
INTEGER*2 IVAR( 3000), IXX( 3000)
C
CHARACTER*70 STR0
CHARACTER*6 STR1( 3)
C
DATA STR0 /'*** F2D-FEM *** ( Max. rel. var. of Velocity (%))'/
DATA STR1 /'0 ','LOP/2',' LOP'/
C
IF ( KEY.EQ.2 ) GO TO 111
C
ICR = 15
CALL SET_TEXT_ATTRIBUTE@( 3, 1., 0., 0. )
CALL DRAW_TEXT@( STR0, 70, 20, ICR )
CALL DRAW_TEXT@( STR1( 2), 290, 450, ICR )
CALL DRAW_TEXT@( STR1( 3), 580, 450, ICR )
C
IX0 = 50
IY0 = 40
CALL RECTANGLE@( IX0, IY0, IX0+580, IY0+400, ICR )
C
111 CONTINUE
IF ( KEY.EQ.1) ICR = 15
IF ( KEY.EQ.2) ICR = 14
IF ( KEY.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 )
CALL DRAW_TEXT@( '(%)', 15, 245, ICR )
C
133 CONTINUE
RIWD = 580.D0/ DFLOAT( LOP)
C
C --------------------------------------
DO 100 I = 1, LOP
C
VRV( I) = DLOG10( VRV( I))
C ********
C
100 CONTINUE
C --------------------------------------
C
SCL = 80.D0
C ------------------------------------------------
DO 200 I = 1, LOP
C
IXX ( I) = IX0 + IDNINT( DFLOAT( I)* RIWD)
IVAR( I) = 280 - IDNINT( VRV( I)* SCL)
C
200 CONTINUE
C ------------------------------------------------
C
*** WRITE(6,*) ' * IVAR( LOP)=', IVAR( LOP)
C
C ------------------------------------------------
C
CALL POLYLINE@( IXX, IVAR, LOP, ICR )
C
C ------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE CNTF ( X, Y, NDE, PV, NNP, NEL )
C
C ----- Contour Lines for Rectangular Meshes ( Linear Interpolation ) --
C
IMPLICIT REAL*8 (A-H,O-Z)
C
INTEGER*2 ICR, IX1, IY1, IX2, IY2
C
COMMON /CNTL/ CON
C
C ----- Important ---
C
DIMENSION X ( NNP), Y ( NNP), PV ( NNP), NDE( NEL,4)
DIMENSION PX( 4), PY( 4), HVN( 10000)
C
ICR = 15
C
IWD = 160
IX0 = 460
IY0 = 70
IX1 = IX0 + IWD
IY1 = IY0 + IWD
CALL RECTANGLE@( IX0, IY0, IX1, IY1, ICR )
C
X0 = 460.D0
Y0 = 70.D0
ICR = 14
SCL = 160.D0
NM = 101
C
TMIN = 100.D0
TMAX = -100.D0
C
C -----------------------------------------------
DO 100 I = 1, NNP
C
IF ( PV( I).GE.TMAX) TMAX = PV( I)
IF ( PV( I).LE.TMIN) TMIN = PV( I)
C
100 CONTINUE
C -----------------------------------------------
C
DTV = 0.02D0 ! For sample running
DTV = 0.05D0 ! For sample running
C
C ------------------------------------------------
KSTR = NM
C
DO 200 I = 1, NM
C
HVN( I) = TMIN + DTV* DFLOAT( I-1)
IF ( HVN( I).GT.TMAX) KSTR = I
IF ( HVN( I).GT.TMAX) GO TO 150
C
200 CONTINUE
C ------------------------------------------------
C
150 CONTINUE
NM = KSTR
C
** WRITE(6,2000) TMAX, TMIN, DTV, NM
** 2000 FORMAT(/,' * Tmax.=',F7.3,' Tmin.=',F7.3,' DTV =',F6.2,
** 1 ' NM =', I3 )
** WRITE(6,*) ' * HVN( 1)=', HVN( 1),' * HVN( NM)=', HVN( NM)
C
C ----- Plot Contour Lines ---
C
PXV = 0.D0
PYV = 0.D0
C
C --------------------------------------------------------------------
DO 300 L = 1, NEL
C
C ----------------------------------------------------------
DO 400 IH = 1, NM
C
HV = HVN( IH)
IC = 0
C
C ------------------------------------------------
DO 500 J = 1, 4
C
I1 = NDE( L,J)
IF ( J .NE.4) THEN
I2 = NDE( L,J+1)
ELSE
I2 = NDE( L,1)
END IF
C
IF ( DABS( PV( I1) - PV( I2)).LE.0.0001D0 ) THEN
IF ( DABS( PV( I1) - HV) .LE.0.0001D0 ) THEN
C
IC = IC + 1
PX( IC) = X( I1)* SCL
PY( IC) = Y( I1)* SCL
END IF
C
ELSE IF ( ( ( PV( I1).LE.HV ) .AND. ( HV.LT.PV( I2)) )
1 .OR.( ( PV( I1).GE.HV ) .AND. ( HV.GT.PV( I2)) ))
2 THEN
C
IC = IC + 1
RT = ( HV - PV( I1))/ ( PV( I2) - PV( I1))
PX( IC) = ( X( I1) + ( X( I2) - X( I1))* RT )* SCL
PY( IC) = ( Y( I1) + ( Y( I2) - Y( I1))* RT )* SCL
END IF
C
500 CONTINUE
C -----------------------------------------------
C
IF ( IC.GE.2) THEN
C
C -----------------------------------------------
DO 600 J = 1, IC
C
JM = J
D2 = ( PXV - PX( J))** 2 + ( PYV - PY( J))** 2
C
C --------------------------------------
DO 650 JJ = J + 1, IC
C
D2V = ( PXV - PX( JJ))** 2 + ( PYV - PY( JJ))** 2
IF ( D2V .LT. D2 ) THEN
JM = JJ
D2 = D2V
END IF
C
650 CONTINUE
C --------------------------------------
PXV = PX( JM)
PYV = PY( JM)
PX( JM) = PX( J)
PY( JM) = PY( J)
PX( J) = PXV
PY( J) = PYV
C
600 CONTINUE
C ------------------------------------------------
C
PXV = PX( IC)
PYV = PY( IC)
C
END IF
C
C ----------------------------------------------------------
IF ( IC.GT.0 ) THEN
C ------------------------------------------------
DO 700 J = 2, IC
C
IX1 = NINT( X0 + PX( J-1))
IY1 = NINT( Y0 + SCL - PY( J-1))
IX2 = NINT( X0 + PX( J))
IY2 = NINT( Y0 + SCL-PY( J))
CALL DRAW_LINE@( IX1, IY1, IX2, IY2, ICR )
C
700 CONTINUE
C ------------------------------------------------
IX1 = NINT( X0 + PX( J))
IY1 = NINT( Y0 - PY( J))
IX2 = NINT( X0 + PX( 1))
IY2 = NINT( Y0 - PY( 1))
IF ( IC.EQ.4 ) CALL DRAW_LINE@( IX1, IY1, IX2, IY2, ICR )
C
END IF
C ----------------------------------------------------------
400 CONTINUE
C ----------------------------------------------------------
C
300 CONTINUE
C --------------------------------------------------------------------
C
RETURN
END
C **********************************************************************