戻る

C
PROGRAM MAIN
C
C ********************************************************************* C
C C
C Graphics for HF2D - FEM C
C C
C ( 2D-Thermal & Fluid Flow Analysis ) C
C C
C ( For Quardrilateral Element ) C
C C
C ( Version 5.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
C
C --- Attention NNP = < 2601, NEL = < 2500 -- Str. Meshes = < 50 X 50 --
C
DIMENSION X( 2601), Y ( 2601), UN( 2601), VN( 2601),
1 T( 2601), NDE( 2500,4)
C ----------------------------------------------------------------------
C
COMMON /CST/ DT, LPM, LOP1, LOP, ALP, NPT, RA, NNP, NEL
C
C ----- Attention --- For Plotting Routines = < 30000 ---
C
COMMON /VRP/ VRV( 30000), VNA( 30000), VNB( 30000), VDN( 30000)
C
C ----- Unstructured Meshes --------------------------------------------
C
** OPEN ( 7, FILE='HF2D_FEM_UE30_RA10_5_R' ) ! DT = 0.0001
** OPEN ( 8, FILE='HF2D_FEM_UE30_RA10_5_Nu')
C
C ----- Structured Meshes ----------------------------------------------
C
** OPEN ( 7, FILE='HF2D_FEM_M50_RA710_R' ) ! DT = 0.001
** OPEN ( 8, FILE='HF2D_FEM_M50_RA710_Nu' )
C
** OPEN ( 7, FILE='HF2D_FEM_M50_RA7100_R' ) ! DT = 0.001
** OPEN ( 8, FILE='HF2D_FEM_M50_RA7100_Nu')
C
OPEN ( 7, FILE='HF2D_FEM_M50_RA10_5_R' ) ! DT = 0.0001
OPEN ( 8, FILE='HF2D_FEM_M50_RA10_5_Nu')
C
** OPEN ( 7, FILE='HF2D_FEM_M50_RA10_6_R' ) ! DT = 0.000 01
** OPEN ( 8, FILE='HF2D_FEM_M50_RA10_6_Nu')
C
** OPEN ( 7, FILE='HF2D_FEM_M50_RA10_7_R' ) ! DT = 0.000 001
** OPEN ( 8, FILE='HF2D_FEM_M50_RA10_7_Nu')
C
C ----------------------------------------------------------------------
C
CALL TTLE
C
CALL PRSB ( X, Y, UN, VN, T, NDE, KEY, ICR )
C
CALL VGA@
C
IF ( KEY.EQ.2) CALL PLNU ( VNA, VNB, 1 )
C
IF ( KEY.EQ.2) GO TO 990
C
CALL PMSH ( X, Y, NDE )
C
CALL PLTN ( X, Y, UN, VN )
C
CALL PLVR ( VRV, 1 )
C
*** CALL PLVR ( VDN, 2 ) ! Attention
C
CALL PCNT ( X, Y, NDE, T )
C
990 CONTINUE
C
CLOSE ( 7)
CLOSE ( 8)
C
STOP
END
C **********************************************************************
C
SUBROUTINE OTLN ( NLIN, ILN, X, Y, NNP )
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
C
DIMENSION ILN( NLIN), X( NNP), Y( NNP)
C
ICR = 15
C
IX0 = 100
IY0 = 70
IWD = 160
C
C ------------------------------------------------
DO 100 I = 1, NLIN
C
N = ILN(( I-1)* 2 + 1)
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 PCNT ( X, Y, NDE, TEP )
C
C ----- To plot contour lines for rectangular meshes ---
C
C ( Linear interpolation ) ---
C
IMPLICIT REAL*8(A-H,O-Z)
C
INTEGER*2 X0, Y0, SCX, SCY, ICR, PX, PY
C
COMMON /CST/ DT, LPM, LOP1, LOP, ALP, NPT, RA, NNP, NEL
C
DIMENSION X ( NNP), Y ( NNP), TEP( NNP), NDE( NEL,4),
1 PX( 4), PY( 4), HVN( 1000)
C
ICR = 15
IWD = 160
C
IX0 = 450
IY0 = 70
IX1 = IX0 + IWD
IY1 = IY0 + IWD
C
CALL RECTANGLE@ ( IX0, IY0, IX1, IY1, ICR )
C
X0 = IX0
Y0 = IY0
ICR = 14
C
SCX = 160
SCY = 160
C
NM = 10
TMN = 0.D0
TMX = 1.D0
C ----------------------------------------------------------
DO 50 I = 1, NM
C
HVN( I) = TMN + ( TMX - TMN)* DFLOAT( I)/ DFLOAT( NM)
IF ( HVN( I).LT.0) NMZ = I
C
50 CONTINUE
C ----------------------------------------------------------
C
C ----- To plot contour lines ---
C
PXV = 0.D0
PYV = 0.D0
C --------------------------------------------------------------------
DO 100 L = 1, NEL
C
C ----------------------------------------------------------
DO 200 IH = 1, NM - 1
C
HV = HVN( IH)
ICL = 0
C -----------------------------------------------------
DO 300 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
C --------------------------------------------------------
IF ( DABS( TEP( I1) - TEP( I2)).LE.0.0001D0) THEN
C
C ---------------------------------------------------
IF ( DABS( TEP( I1) - HV ).LE.0.0001D0) THEN
ICL = ICL + 1
PX( ICL) = X( I1)* SCX
PY( ICL) = Y( I1)* SCY
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))))
2 THEN
C
ICL = ICL + 1
RT = ( HV - TEP( I1))/ ( TEP( I2) - TEP( I1))
PX( ICL) = ( X( I1) + ( X( I2) - X( I1))* RT)* DFLOAT( SCX)
PY( ICL) = ( Y( I1) + ( Y( I2) - Y( I1))* RT)* DFLOAT( SCY)
C
END IF
C ----------------------------------------------------
C
300 CONTINUE
C ----------------------------------------------------
IF ( ICL.GE.2) THEN
C
C ------------------------------------------------
DO 400 J = 1, ICL
C
JM = J
D2 = ( PXV - PX( J))** 2 + ( PYV - PY( J))** 2
C
C --------------------------------------
DO 500 JJ = J+1, ICL
C
D2V = ( PXV - PX( JJ))** 2 + ( PYV - PY( JJ))** 2
C
IF ( D2V.LT.D2) THEN
JM = JJ
D2 = D2V
END IF
C
500 CONTINUE
C --------------------------------------
PXV = PX( JM)
PYV = PY( JM)
PX(JM) = PX( J)
PY(JM) = PY( J)
PX( J) = PXV
PY( J) = PYV
C
400 CONTINUE
C ------------------------------------------------
PXV = PX( ICL)
PYV = PY( ICL)
END IF
C --------------------------------------
C
C ------------------------------------------------
IF ( ICL.GT.0 ) THEN
C ------------------------------------------------
DO 600 J = 2, ICL
C
CALL DRAW_LINE@ ( X0+PX( J-1), Y0+SCY-PY( J-1),
1 X0+PX( J), Y0+SCY-PY( J), ICR )
C
600 CONTINUE
C ------------------------------------------------
IF ( ICL.EQ.4)
1 CALL DRAW_LINE@ ( X0+PX( J), Y0-PY( J),
2 X0+PX( 1), Y0-PY( 1), ICR )
END IF
C ------------------------------------------------
C
200 CONTINUE
C ----------------------------------------------------------
C
100 CONTINUE
C --------------------------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE PLNU ( VNA, VNB, KEY )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
INTEGER*2 ICR, IX0, IY0
C
COMMON /CST/ DT, LPM, LOP1, LOP, ALP, NPT, RA, NNP, NEL
C
C ----- Attention N.LE.30000 ---
C
DIMENSION VNA( LOP), VNB( LOP)
C
INTEGER*2 IVRA( 30000), IVRB( 30000), IXX( 30000)
C
CHARACTER*70 STR0
CHARACTER*6 STR1( 3)
CHARACTER*4 STR2( 4)
C
DATA STR0 /'*** HF2D - FEM (uvp) *** ( Value of Nu )'/
DATA STR1 /'0 ','LOP/2',' LOP'/
DATA STR2 /' 20',' 15',' 10',' 5'/
C
IF ( KEY.EQ.2) GO TO 111
C
ICR = 15
C
CALL SET_TEXT_ATTRIBUTE@ ( 3, 1., 0., 0. )
CALL DRAW_TEXT@ ( STR0, 70, 20, ICR )
C
CALL SET_TEXT_ATTRIBUTE@ ( 3, 1.5, 0., 0. )
C
IF ( DABS( RA-7100.D0 ).LE.0.0001D0 )
1 CALL DRAW_TEXT@ ( ' Ra = 7100', 290, 320, ICR )
C
IF ( DABS( RA-710.D0 ).LE.0.0001D0 )
1 CALL DRAW_TEXT@ ( ' Ra = 710', 290, 320, ICR )
C
IF ( DABS( RA-100000.D0 ).LE.0.0001D0 )
1 CALL DRAW_TEXT@ ( ' Ra = 10**5', 200, 300, ICR )
C
IF ( DABS( RA-1.0D+6 ).LE.0.0001D0 )
1 CALL DRAW_TEXT@ ( ' Ra = 10**6', 70, 300, ICR )
C
IF ( DABS( RA-1.0D+7 ).LE.0.0001D0 )
1 CALL DRAW_TEXT@ ( ' Ra = 10**7', 200, 300, ICR )
C
C ----------------------------
C
CALL SET_TEXT_ATTRIBUTE@ ( 3, 1., 0., 0. )
ICR = 15
C
CALL DRAW_TEXT@ ( STR1( 1), 50, 450, ICR )
CALL DRAW_TEXT@ ( STR1( 2), 315, 450, ICR )
CALL DRAW_TEXT@ ( STR1( 3), 580, 450, ICR )
C
CALL DRAW_TEXT@ ( STR2( 1), 15, 40, ICR )
CALL DRAW_TEXT@ ( STR2( 2), 15, 140, ICR )
CALL DRAW_TEXT@ ( STR2( 3), 15, 240, ICR )
CALL DRAW_TEXT@ ( STR2( 4), 15, 340, ICR )
C
IX0 = 50
IY0 = 50
CALL RECTANGLE@ ( IX0, IY0, IX0+580, IY0+400, ICR )
C
111 CONTINUE
C
IF ( KEY.EQ.1) ICR = 15
RIWD = 580.D0/ DFLOAT( LOP)
C
C ----------------------------------------------------------
DO 100 I = 1, LOP
C
IXX ( I) = IX0 + IDNINT( DFLOAT( I)* RIWD)
IVRA( I) = IY0 + 400 - IDNINT( VNA( I)* 20.D0)
IVRB( I) = IY0 + 400 - IDNINT( VNB( I)* 20.D0)
C
100 CONTINUE
C ----------------------------------------------------------
C
ICR = 12
CALL POLYLINE@ ( IXX, IVRA, LOP, ICR )
C
ICR = 14
CALL POLYLINE@ ( IXX, IVRB, LOP, ICR )
C
RETURN
END
C **********************************************************************
C
SUBROUTINE PLTN ( X, Y, UN, VN )
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 /CST/ DT, LPM, LOP1, LOP, ALP, NPT, RA, NNP, NEL
C
SCX = 240.D0
SCY = 240.D0
C
WRITE(6,2000) RA
2000 FORMAT(' * Ra =',F9.2)
C
IF ( RA.LT.7000.D0) SCL = 2.0D0
IF ( RA.GE.7000.D0) SCL = 0.5D0
IF ( RA.GE.1.0D5) SCL = 0.1D0
IF ( RA.GE.1.0D6) SCL = 0.05D0
C
ICR = 15
IX0 = 270
IY0 = 70
IWD = 160
TL = DT* DFLOAT( LOP)
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 + 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 PLVR ( VRV, KEY )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
INTEGER*2 ICR, IX0, IY0
C
C ------ Attention N.LE.30000 ---
C
DIMENSION VRV( LOP)
C
INTEGER*2 IVAR( 30000), IXX( 30000)
C
COMMON /CST/ DT, LPM, LOP1, LOP, ALP, NPT, RA, NNP, NEL
C
CHARACTER*110 STR0
CHARACTER*6 STR1( 3)
C
DATA STR0 /'** HF2D-FEM (uvp) *(Max.rel.var.of Velocity(%) & Diff
1. of Nu-Value(-))'/
DATA STR1 /'0 ','LOP/2',' LOP'/
C
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 SET_TEXT_ATTRIBUTE@ ( 3, 1.5, 0., 0. )
C
IF ( DABS( RA-7100.D0 ).LE.0.0001D0 )
1 CALL DRAW_TEXT@ ( ' Ra = 7100', 150, 370, ICR )
C
IF ( DABS( RA-710.D0 ).LE.0.0001D0 )
1 CALL DRAW_TEXT@ ( ' Ra = 710', 150, 370, ICR )
C
IF ( DABS( RA-100000.D0 ).LE.0.0001D0 )
1 CALL DRAW_TEXT@ ( ' Ra = 10**5', 130, 370, ICR )
C
IF ( DABS( RA-1.0D+6 ).LE.0.0001D0 )
1 CALL DRAW_TEXT@ ( ' Ra = 10**6', 70, 370, ICR )
C
IF ( DABS( RA-1.0D+7 ).LE.0.0001D0 )
1 CALL DRAW_TEXT@ ( ' Ra = 10**7', 70, 400, ICR )
C
C ----------------------------
C
CALL SET_TEXT_ATTRIBUTE@ ( 3, 1., 0., 0. )
C
ICR = 15
CALL DRAW_TEXT@ ( STR1( 1), 40, 450, ICR )
CALL DRAW_TEXT@ ( STR1( 2), 290, 450, ICR )
CALL DRAW_TEXT@ ( STR1( 3), 580, 450, ICR )
IX0 = 50
IY0 = 40
CALL RECTANGLE@ ( IX0, IY0, IX0+580, IY0+400, ICR )
C
111 CONTINUE
C
IF ( KEY.NE.1) GO TO 133
C
CALL DRAW_TEXT@ ( '1000', 15, 38, ICR )
CALL DRAW_TEXT@ ( ' 100', 15, 118, ICR )
CALL DRAW_TEXT@ ( ' 10', 15, 198, ICR )
CALL DRAW_TEXT@ ( ' 1', 15, 278, ICR )
CALL DRAW_TEXT@ ( ' 0.1', 15, 358, ICR )
CALL DRAW_TEXT@ ( '0.01', 15, 438, ICR )
CALL DRAW_TEXT@ ( '(%)' , 15, 245, ICR )
C
133 CONTINUE
IX0 = 50
IY0 = 40
CALL RECTANGLE@ ( IX0, IY0, IX0 + 580, IY0 + 400, ICR )
RIWD = 580.D0/ DFLOAT( LOP)
C
C ------------------------------------------------
DO 100 I = 1, LOP
C
IF ( VRV( I).LE.0.01D0 ) GO TO 333 ! Attention
C *********
C
VRV( I) = DLOG10( VRV( I))
C
100 CONTINUE
333 LOP = I - 1
WRITE(6,1233) LOP
1233 FORMAT(' *** Attention : Revised No. of Plottings =', I8)
C ------------------------------------------------
C
SCL = 80.D0
C ------------------------------------------------
IF ( KEY.EQ.2) ICR = 14
C
RIWD = 580.D0/ DFLOAT( LOP)
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
CALL POLYLINE@ ( IXX, IVAR, LOP, ICR )
C
C ---------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE PMSH ( X, Y, NDE )
C
IMPLICIT REAL*8(A-H,O-Z)
C
INTEGER*2 IX1, IX2, IY1, IY2, ICR
C
COMMON /CST/ DT, LPM, LOP1, LOP, ALP, NPT, RA, NNP, NEL
C
DIMENSION X( NNP), Y( NNP), NDE( NEL,4)
C
CX = 90.D0
CY = 230.D0
C
C ------------------------------------------------
DO 100 I = 1, NEL
C
C ------------------------------------------------
DO 100 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
100 CONTINUE
C ------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE PRSB ( X, Y, UN, VN, T, NDE, KEY, ICR )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
INTEGER*2 ICR, ICN
C
DIMENSION X ( NNP), Y( NNP), UN( NNP), VN( NNP), T( NNP),
1 NDE( NEL,4)
C
COMMON /CST/ DT, LPM, LOP1, LOP, ALP, NPT, RA, NNP, NEL
COMMON /VRP/ VRV( 30000), VNA( 30000), VNB( 30000), VDN( 30000)
C
ICR = 15
ICN = 1
C
WRITE(6,2000) RA, DT, NNP, NEL, NPT
2000 FORMAT(' * Ra =', F10.1,' * DT =',F7.5,' NNP =', I6,
1 ' NEL =', I6,' NPT =', I3 )
C
READ(7,1000) ( X( I), Y( I), I=1, NNP)
1000 FORMAT(10F8.4)
READ(7,1100) (( NDE( I,J), J=1, 4), I= 1, NEL)
1100 FORMAT(20I4)
READ(7,1200) LOP
1200 FORMAT(I10)
WRITE(6,*) LOP
C
C ----------------------------------------------------------
XMN = 1.D0
C
DO 100 I = 2, NNP
C
XX = X( I)
IF ( DABS( XX).GT.0.0001D0) XMN = MIN( XMN, XX)
100 CONTINUE
C ----------------------------------------------------------
WRITE(6,2100) LOP, XMN
2100 FORMAT(' * Loop from HF2D-FEM =',I6,' XMN.=',F7.4)
C
C ----- Last Loop ---
C
READ(7,1300) ( VRV( I), I = 1, LOP)
1300 FORMAT(10F11.4)
*** WRITE(6,*) ( VRV( I), I = 1, LOP)
C
READ(7,1300) ( VDN( I), I = 1, LOP)
READ(7,1300) ( UN( I), VN( I), I = 1, NNP)
READ(7,1300) ( T( I), I = 1, NNP )
READ(8,1000) ( VNA( I), VNB( I), I = 1, LOP)
C ----------------------------------------------------------
C
C ----- Attention --- Adjusting for No. of Plottings ---
C ===> 1000 plottings ---
C
C ------------------------------------------------
IF ( LOP.GE.30000) THEN
C
NL = 0
C --------------------------------------
DO 200 I = 1, LOP, 30
C
NL = NL + 1
VRV( NL) = VRV( I)
VDN( NL) = VDN( I)
VNA( NL) = VNA( I)
VNB( NL) = VNB( I)
C
200 CONTINUE
C --------------------------------------
WRITE(6,*) ' *** Attention: No_Plot_LOP =', NL
LOP = NL
C
END IF
C ------------------------------------------------
WRITE(6,*) ' '
WRITE(6,*) ' * Key = 1 : Meshes, Vel.-Distr., Temp-Distr. & Var.of
1 Vel. etc.'
WRITE(6,*) ' = 2 : Plot of Nu-Value '
READ (5,*) KEY
WRITE(6,*) ' *** KEY =', KEY
C
RETURN
END
C **********************************************************************
C
SUBROUTINE TTLE
C
IMPLICIT REAL*8(A-H,O-Z)
C
COMMON /CST/ DT, LPM, LOP1, LOP, ALP, NPT, RA, NNP, NEL
C
WRITE(6,*) ' '
WRITE(6,*) 'C ************************************************* C'
WRITE(6,*) 'C C'
WRITE(6,*) 'C Graphics for HF2D - FEM (uvp) C'
WRITE(6,*) 'C C'
WRITE(6,*) 'C ( Version 5.0 ) C'
WRITE(6,*) 'C C'
WRITE(6,*) 'C 2D Thermal Fluid Analysis by FEM C'
WRITE(6,*) 'C C'
WRITE(6,*) 'C ( For Quardrilateral 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
READ(7,1000) DT, RA, NNP, NEL, NPT
1000 FORMAT(F12.9, F12.1, 3I5)
C
RETURN
END
C **********************************************************************