戻る

C
PROGRAM MAIN
C *******************************************************************C
C C
C HF2D - FET( SV ) - FIG C
C C
C ( V.4.3 ) C
C C
C ( Thermal Fluid Analysis Software by Finite Element Method ) C
C C
C ( Triangular-Linear Element ) C
C C
C Copyright : Yasuhiro MATSUDA C
C C
C *******************************************************************C
C
IMPLICIT REAL*8(A-H,O-Z)
C
INTEGER*2 ICR, NRS, ICN
C
C ----- Attention : DATA for ( 17 x 17 x 4 : 113 / 196 ) ---
C ==================================================
C
DIMENSION X ( 113), Y ( 113), UN ( 113), VN( 113),
1 PSI( 113), VOR( 113), TEP( 113)
C
INTEGER*2 NDE, IRS
DIMENSION NDE( 196,3), IRS( 616)
C
COMMON /CNST/ DT, LPLT, LOP1, LOP, ALP
C
C ----- Attention : < = 2000 ---
COMMON /VRPT/ VARV( 2000), VRT( 2000), VRP( 2000)
C
OPEN (10, FILE='HFE_T_FIG.DAT' )
C
CALL TTLE
C
READ(10,1000) ICN, LOP, DT, RE, RA, PRN, NNP, NEL
1000 FORMAT(I2, I4, F6.3, 3F8.3, 2I4)
C
CALL PREP ( PSI, VOR, TEP, X, Y, UN, VN, NDE, IRS,
1 NRS, ICN, ICR, NNP, NEL, ICN )
C
CALL VGA@
C
CALL OTLN ( NRS, IRS, X, Y, NNP, ICN )
C
CALL PLTN ( X, Y, UN, VN, NNP, ICN )
C
CALL PLVR ( VARV, LOP, 1, ICN )
C
IF ( ICN.EQ.0) CALL PLVR ( VRP, LOP, 3, ICN )
C
IF ( ICN.EQ.1) THEN
C
CALL PLVR ( VRT, LOP, 4, ICN )
C
END IF
C
CALL PLCNT ( X, Y, NDE, PSI, NNP, NEL, 0 )
C
IF ( ICN.EQ.1) CALL PLCNT ( X, Y, NDE, TEP, NNP, NEL, 1 )
C
CLOSE ( 10)
C
STOP
END
C **********************************************************************
C
SUBROUTINE PREP ( PSI, VOR, TEP, X, Y, UN, VN, NDE,
1 IRS, NRS, ICN, ICR, NNP, NEL )
C
IMPLICIT REAL*8(A-H,O-Z)
C
INTEGER*2 ICR, NRS, ICN, NDE, IRS
C
COMMON /CNST/ DT, LPLT, LOP1, LOP, ALP
COMMON /VRPT/ VARV( 2000), VRT( 2000), VRP( 2000)
C
DIMENSION PSI( NNP), VOR( NNP), TEP( NNP), X ( NNP),
1 Y ( NNP), UN ( NNP), VN ( NNP), NDE( NEL,3),
2 IRS( 616)
C
REWIND 10
C
READ(10,1000) ICN, LOP, DT, RE, RA, PRN
1000 FORMAT(I2, I4, F6.3, 3F8.1 )
C
IF ( ICN.EQ.1)
1WRITE(6,2000) ICN, LOP, DT, RA, PRN, NNP, NEL
2000 FORMAT(' * ICN =',I2,' LOP =',I3,' DT =',F8.4,' Ra =',F9.1,
1 ' Pr =',F5.1,/,12X,' NNP =',I4,' NEL =',I4,/)
IF ( ICN.EQ.0)
1WRITE(6,2100) ICN, LOP, DT, RE, NNP, NEL
2100 FORMAT(' * ICN =',I2,' LOP =',I3,' DT =',F8.4,' Re =',F6.1,
1 /,10X,' NNP =',I4,' NEL =',I4,/)
C
ICR = 15
C
READ(10,1100) ( X( I), Y( I), I = 1, NNP)
1100 FORMAT(10F8.4)
READ(10,1200) ( ( NDE( I,J ), I = 1, NEL ), J = 1, 3)
1200 FORMAT(20I4)
READ(10,1200) NRS, ( IRS( I), I=1, 616 )
C
C ----- At the Last LOP ---
C
READ(10,1300) ( VRP( I), I = 1, LOP )
1300 FORMAT(10F12.5)
C
IF ( ICN.EQ.1) READ(10,1300) ( VARV( I), VRT( I), I = 1, LOP )
IF ( ICN.EQ.0) READ(10,1300) ( VARV( I), I = 1, LOP )
READ(10,1100) ( UN( I), VN( I), I = 1, NNP )
C
IF ( ICN.EQ.1)
1 READ(10,1100) ( PSI( I), VOR( I), TEP( I), I = 1, NNP )
IF ( ICN.EQ.0) READ(10,1100) ( PSI( I), VOR( I), I = 1, NNP )
C
RETURN
END
C **********************************************************************
C
SUBROUTINE TTLE
C
IMPLICIT REAL*8(A-H,O-Z)
C
WRITE(6,*) ' '
WRITE(6,*) 'C ************************************************* C'
WRITE(6,*) 'C C'
WRITE(6,*) 'C H F 2 D - F E T ( SV ) - FIG C'
WRITE(6,*) 'C C'
WRITE(6,*) 'C ( Version 4.3 ) C'
WRITE(6,*) 'C C'
WRITE(6,*) 'C 2D Thermal Fluid Analysis by FEM C'
WRITE(6,*) 'C C'
WRITE(6,*) 'C ( Triangular Linear Element ) C'
WRITE(6,*) 'C C'
WRITE(6,*) 'C Copyright 2014 : Yasuhiro MATSUDA C'
WRITE(6,*) 'C C'
WRITE(6,*) 'C ************************************************* C'
WRITE(6,*) ' '
C
RETURN
END
C **********************************************************************
C
SUBROUTINE PLCNT ( X, Y, NDE, TEP, NNP, NEL, KKF )
C
C ----- To plot CONTour lines for triangular meshes
C ( Linear interpolation ) ---
C
IMPLICIT REAL*8(A-H,O-Z)
C
INTEGER*2 NDE
INTEGER*2 ICR, IWD, IX0, IY0, IX1, IY1, IX2, IY2, IX3, IY3
C
DIMENSION X ( NNP), Y ( NNP), NDE( NEL,3), TEP( NNP), PX( 3),
1 PY( 3), HVN( 50)
C
IWD = 160
C
IF ( KKF.EQ.0) IX0 = 100
IF ( KKF.EQ.0) ICR = 11
C
IF ( KKF.EQ.1) IX0 = 460
IF ( KKF.EQ.1) ICR = 12
C
IY0 = 70
IX1 = IX0 + IWD
IY1 = IY0 + IWD
C
CALL RECTANGLE@( IX0, IY0, IX1, IY1, 15 )
C
NM = 10
C
C ---------------------------------------------------------------
C
TMIN = 100.D0
TMAX = -100.D0
C ------------------------------------------------
DO 100 I = 1, NNP
C
IF ( TEP( I).GE.TMAX) TMAX = TEP( I)
IF ( TEP( I).LE.TMIN) TMIN = TEP( I)
C
100 CONTINUE
C ------------------------------------------------
WRITE(6,2000) TMAX, TMIN
2000 FORMAT(' * Cmax.=',F7.3,' * Cmin.=',F7.3)
C
C ------------------------------------------------
RWD = ( TMAX - TMIN )/ DFLOAT( NM)
C
C ------------------------------------------------
DO 200 I = 1, NM
C
HVN( I) = TMIN + RWD* DFLOAT( I)
IF ( HVN( I).LT.0.D0 ) NMZ = I
C
200 CONTINUE
C ------------------------------------------------
C
C ----- To Plot Contour Lines ---
C
PXV = 0.D0
PYV = 0.D0
C
C ----------------------------------------------------------
DO 300 L = 1, NEL
C
C -----------------------------------------------
DO 400 IH = 1, NM-1
C
HV = HVN( IH)
ICL = 0
C
C --------------------------------------
DO 500 J = 1, 3
C
I1 = NDE( L,J )
IF ( J .NE. 3) THEN
I2 = NDE( L,J+1 )
ELSE
I2 = NDE( L,1 )
END IF
C
C ----------------------------------------------------------
IF ( DABS ( TEP(I1)-TEP(I2)) .LE. 0.0001D0 ) THEN
C ------------------
IF ( DABS ( TEP(I1)-HV ) .LE. 0.0001D0 ) THEN
ICL = ICL + 1
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))
PX(ICL) = ( X(I1) + ( X(I2) - X(I1))*RT ) * DFLOAT( IWD)
PY(ICL) = ( Y(I1) + ( Y(I2) - Y(I1))*RT ) * DFLOAT( IWD)
END IF
C ----------------------------------------------------------
C
500 CONTINUE
C --------------------------------------
C
IF (ICL .GE. 2) THEN
C
C --------------------------------------
DO 600 J = 1, ICL
C
JM = J
D2 = ( PXV - PX( J))** 2 + ( PYV - PY( J))** 2
C
C ----------------------------
DO 700 JJ = J+1, ICL
C
D2V = ( PXV - PX( JJ))** 2 + ( PYV - PY( JJ))** 2
IF ( D2V .LT. D2 ) THEN
JM = JJ
D2 = D2V
END IF
C
700 CONTINUE
C ----------------------------
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(ICL)
PYV = PY(ICL)
END IF
C
C --------------------------------------
IF ( ICL.GT.0 ) THEN
C
C ------------------------------------------------
DO 800 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
800 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 ----------------------------
C
END IF
C --------------------------------------
C
400 CONTINUE
C ------------------------------------------------
C
300 CONTINUE
C ----------------------------------------------------------
C
RETURN
END
C *******************************************************************
C
SUBROUTINE PLTN ( X, Y, UN, VN, NNP, ICN )
C
IMPLICIT REAL*8(A-H,O-Z)
C
INTEGER*2 ICN, ICR, IWD, IX0, IY0, IX1, IY1, IX2, IY2, GX, GY
C
DIMENSION X( NNP), Y( NNP), UN( NNP), VN( NNP)
C
COMMON /CNST/ DT, LPLT, LOP1, LOP, ALP
C
SCLX = 240.D0
SCLY = 240.D0
SCL = 50.D0
C
IF ( ICN.EQ.1 ) SCL = 10.D0
ICR = 15
C
IX0 = 280
IY0 = 70
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
C ----------------------------------------------------------
DO 100 I = 1, NNP
C
GX = IX0 + IDNINT( X( I)* DFLOAT( IWD))
GY = IY0 + IWD - IDNINT( 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 ----------------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE OTLN ( NLIN, ILN, X, Y, NNP, ICN )
C
C ----- Fig. of Meshes ---
C
IMPLICIT REAL*8(A-H,O-Z)
C
INTEGER*2 ILN, NLIN, ICR, IWD, IX0, IY0, IX1, IY1, IX2, IY2, ICN
C
DIMENSION ILN( NLIN), X( NNP), Y( NNP)
C
ICR = 15
C
IX0 = 460
IY0 = 250
IWD = 160
C
IF ( ICN.EQ.0) IY0 = 70
C
C ----------------------------------------------------------
DO 100 I = 1, NLIN
C
N = ILN( (I-1)* 2 + 1 )
C
IX1 = IDNINT( X( N)* DFLOAT( IWD)) + IX0
IY1 = IDNINT( Y( N)* DFLOAT( IWD)) + IY0
C
N = ILN( I* 2)
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 ----------------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE PLVR ( RMX, LOP, KEY, ICN )
C
IMPLICIT REAL*8(A-H,O-Z)
C
INTEGER*2 ICR, IX0, IY0
INTEGER*2 ICN
C
C ----- Attention N.LE.500 ---
C
DIMENSION RMX( LOP)

INTEGER*2 IVR( 2000), IXX( 2000)
C
CHARACTER*70 STR0
CHARACTER*70 STR3
CHARACTER*5 STR1(3)
C
DATA STR0 /'*** HF2D-FET (SV) *** ( Max. rel. var. of Vorticity &
1 Temperature (%))'/
DATA STR3 /'*** HF2D-FET (SV) *** ( Max. rel. var. of Str. Functi
1on & Vorticity (%))'/

DATA STR1/'0 ','NPL/2',' NPL'/
C
CALL SET_TEXT_ATTRIBUTE@( 3, 1., 0., 0. )
C
ICR = 15
C
IF ( ICN.EQ.1 ) CALL DRAW_TEXT@( STR0, 70, 20, ICR )
IF ( ICN.EQ.0 ) CALL DRAW_TEXT@( STR3, 70, 20, 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 ( KEY.GE.2 ) GO TO 1234
C
ICR = 15
C
CALL SET_TEXT_ATTRIBUTE@( 3, 1., 0., 0. )
C
IF ( ICN.EQ.1) CALL DRAW_TEXT@( STR0, 70, 20, ICR )
IF ( ICN.EQ.0) CALL DRAW_TEXT@( STR3, 70, 20, 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 ( 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 )
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
RIWD = 580.D0/ DFLOAT( LOP)
C ---------------------------------
C
IF ( KEY.EQ.1 ) ICR = 11
IF ( KEY.EQ.2 ) ICR = 15
IF ( KEY.EQ.3 ) ICR = 9
IF ( KEY.EQ.4 ) ICR = 12
C
C ------------------------------------------------
DO 100 I = 1, LOP
C
IF ( RMX( I).LE.0.D0 ) GO TO 100
RMX( I) = DLOG10( RMX( I))
C
100 CONTINUE
C ------------------------------------------------
DO 200 I = 1, LOP
C
IXX( I) = IX0 + IDNINT( DFLOAT( I)* RIWD )
IVR( I) = 280 - IDNINT( 80.D0* RMX( I))
C
200 CONTINUE
C
C ------------------------------------------------
C
CALL POLYLINE@( IXX, IVR, LOP, ICR )
C
C ------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE PLUV ( UT, VT )
C
IMPLICIT REAL*8(A-H,O-Z)
C
INTEGER*2 XINT, YINT, X0, Y0, X1, Y1, IWD
C
C ------ Attention ---
DIMENSION UT( 24,24), VT( 24,24)
C
COMMON /ABC/ XINT, YINT, IWD
C
SCL = 20.D0
XINT = XINT + IWD + 20
C
CALL RECTANGLE@( XINT, YINT, XINT+IWD, YINT+IWD, 15 )
C
IP = 20
JP = 20
C
IP1 = 23
JP1 = 23
C
XXL = IWD / 20
YYL = IWD / 20
C
C ----------------------------------------------------------
DO 100 I = 1, IP
DO 100 J = 1, JP
C ----------------------------------------------------------
C
C ----- To change of Velocity direction : Conversion of X-Y ---
C
X0 = IDNINT( YYL* ( DFLOAT( I)-0.5D0 )) + XINT
Y0 = IDNINT( XXL* ( DFLOAT( J)-0.5D0 )) + YINT
C
X1 = X0 - IDNINT( ( VT( JP1-I,J+2 )* SCL))
Y1 = Y0 + IDNINT( ( UT( IP1-I,J+2 )* SCL))
C
C- ---- To omit 'PLOT' in case of 'UV'.lt.0.01 ---
C
CALL DRAW_LINE@( X0, Y0, X1, Y1, 15 )
C
100 CONTINUE
C ----------------------------------------------------------
C
RETURN
END
************************************************************************