! copyright 2005, J. E. Akin, all rights reserved. ! file: printing_lib.f SUBROUTINE AT (LINE) ! * * * * * * * * * * * * * * * * * * * * * ! DEBUGGING AID, PRINT AN INTEGER ! * * * * * * * * * * * * * * * * * * * * * INTEGER, INTENT(IN) :: LINE WRITE (6, *) ' -------->>> HERE AT ', LINE END SUBROUTINE AT SUBROUTINE NODAL_BAR_CHART (IBAR, I_PARM, X, D) ! * * * * * * * * * * * * * * * * * * * * * * * * * * ! PRINT-PLOT BAR CHARTS OF NODAL PARAMETER I_PARM AT ! THE NODES IN ARRAY NODE AND SCALE THE RELATIVE ! DISTANCE BETWEEN THE POINTS ! * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants IMPLICIT NONE INTEGER, PARAMETER :: MID = 25, LINE = 2 * MID REAL(DP), INTENT(IN) :: X (MAX_NP, N_SPACE), D (N_D_FRE) INTEGER, INTENT(IN) :: IBAR, I_PARM ! Automatic Arrays INTEGER :: NODE_BAR (IBAR) ! Locals INTEGER :: I, IS, INDEX, IPMNG, ISPACE, J, JY, JZ, & K, KOUNT, L, LIMIT, NLAST, NOPLOT REAL(DP) :: DIST, DMIN, DMINSQ, DTEST, RANGE, YMAX, YMIN CHARACTER ALINE (LINE+1), SKIP (LINE+1) CHARACTER(1), SAVE :: BLANK = ' ', DOT = '+', DASH = '-', & A_X = 'X', PLUS = '+', A_0 = '0' DATA NOPLOT / 0 / ! MAX_NP = TOTAL NUMBER OF NODES IN SYSTEM ! I_PARM = NODAL PARAMETER TO BE GRAPHED,1<=I_PARM<=N_G_DOF ! NODE_BAR = LIST OF NODES TO BE USED (WHEN IBAR > 1) ! NO_DIST = 1 OMIT DISTANCES BETWEEN NODAL BAR LINES ! X = ARRAY OF GLOBAL COORDINATES OF ALL NODES ! N_D_FRE = TOTAL NO. OF DEGREES OF FREEDOM IN SYSTEM ! D = ARRAY OF ALL NODAL PARAMETERS IN THE SYSTEM ! N_G_DOF = NUMBER OF PARAMETERS PER NODE ! N_SPACE = DIMENSION OF SOLUTION SPACE ! IBAR = NUMBER OF NODES TO BE INCLUDED IN BAR CHART ! IF IBAR=1 USE ALL NODES, NODE_BAR NOT USED NOPLOT = NOPLOT + 1 LIMIT = IBAR IF ( IBAR == 1) LIMIT = MAX_NP WRITE (N_PRT, 5000) NOPLOT, DOF_NAME (I_PARM), LIMIT 5000 FORMAT ('*** PRINT PLOT NUMBER',I3,' ***',/, & A12, ' EVALUATED AT ',I5,' NODE POINTS') !b IPMNG = I_PARM - N_G_DOF IF ( NO_DIST == 0) THEN ! FIND MINIMUM DISTANCE BETWEEN POINTS DMINSQ = 0.d0 DO IS = 1, N_SPACE DMINSQ = DMINSQ + (X (2, IS) - X (1, IS) ) **2 END DO DO J = 2, LIMIT I = J IF ( IBAR > 1) I = NODE_BAR (J) DTEST = 0.0 DO IS = 1, N_SPACE DTEST = DTEST + (X (I, IS) - X (I - 1, IS) ) **2 END DO IF ( DTEST < DMINSQ) DMINSQ = DTEST END DO DMIN = SQRT (DMINSQ) END IF !--> ESTABLISH GRAPH LIMITING VALUES ! INITALIZE MAX MIN VALUES OF PARAMETER L = 1 IF ( IBAR /= 1) L = NODE_BAR (1) !b INDEX = N_G_DOF * L + IPMNG INDEX = N_G_DOF * (L - 1) + I_PARM YMAX = D (INDEX) YMIN = YMAX DO I = 1, LIMIT L = I IF ( IBAR /= 1) L = NODE_BAR (I) !b INDEX = N_G_DOF * L + IPMNG INDEX = N_G_DOF * (L - 1) + I_PARM DTEST = D (INDEX) IF ( DTEST > YMAX) YMAX = DTEST IF ( DTEST < YMIN) YMIN = DTEST END DO KOUNT = 1 WRITE (N_PRT, 5010) YMIN, YMAX 5010 FORMAT (' RANGE ON GRAPH IS ',1PE12.5,' TO ',1PE12.5/) IF ( YMIN == YMAX ) THEN WRITE (N_PRT, 5011) 5011 FORMAT ('CONSTANT VALUE, PLOT SKIPPED') RETURN END IF RANGE = YMAX - YMIN ! FLAG ZERO VALUE IF ( YMIN < 0.0 .AND. YMAX > 0.0 ) THEN JZ = 1*(YMAX - 0.0) / RANGE & + (LINE+1)*(0.0 - YMIN) / RANGE + 0.1 ELSE JZ = 0 END IF ! SCALING COMPLETE SKIP (2:LINE) = BLANK ALINE (2:LINE) = DASH SKIP (1) = PLUS SKIP (MID+1) = DOT SKIP (LINE+1) = PLUS ALINE (1:LINE+1:5) = PLUS ALINE (MID+1) = DOT WRITE (N_PRT, 5020) ALINE 5020 FORMAT(' NODE VALUE ', (101A1) ) L = 1 IF ( IBAR > 1) L = NODE_BAR (1) DO K = 1, LIMIT NLAST = L L = K IF ( IBAR > 1) L = NODE_BAR (K) !b INDEX = N_G_DOF * L + IPMNG INDEX = N_G_DOF * (L - 1) + I_PARM DTEST = D (INDEX) ALINE (2:LINE) = BLANK ALINE (1) = PLUS ALINE (LINE+1) = PLUS IF ( JZ > 0) ALINE (JZ) = A_0 ISPACE = 1 IF ( NO_DIST == 0) THEN !--> FIND DISTANCE BETWEEN TWO POINTS DIST = SUM ((X (L, :) - X (NLAST, :))**2) !9 ! MINIMUM DISTANCE IS ONE SPACE DIST = SQRT (DIST) IF ( K /= 1 ) THEN IF ( DMIN > 0.0 ) THEN ISPACE = DIST / DMIN + 0.5 ELSE ISPACE = 1 END IF END IF IF ( ISPACE > 10) ISPACE = 5 END IF !b IF ( NO_DIST > 0 .AND. K > 1 ) THEN IF ( NO_DIST == 0 .AND. K > 1 ) THEN DO I = 1, ISPACE WRITE (N_PRT, 5030) SKIP 5030 FORMAT (20X, (101A1) ) END DO END IF ! LINEAR INTERPOLATION FOR COLUMN NUMBER JY = 1 * (YMAX - DTEST) / RANGE & + (LINE+1) * (DTEST - YMIN) / RANGE + 0.1 ALINE (1:JY) = A_X WRITE (N_PRT, 5040) L, D (INDEX), ALINE 5040 FORMAT (I7,2X,1PE10.3,1X,101A1) END DO ALINE (2:LINE) = DASH ALINE (1:LINE+1:5) = PLUS ALINE (MID+1) = DOT WRITE (N_PRT, 5020) ALINE END SUBROUTINE NODAL_BAR_CHART SUBROUTINE G_PRINT (A, NR, NC) ! ------------------------------------------------------ ! F90 PRINTING OF REAL MATRIX A(NR,NC) WITH G FORMAT ! ------------------------------------------------------ Use Precision_Module IMPLICIT NONE INTEGER, INTENT(IN) :: NC, NR REAL(DP), INTENT(IN) :: A(NR*NC) INTEGER, PARAMETER :: MAX = 7 INTEGER :: NCOL(MAX), I, IC, J, K, L, MAXCOL INTEGER :: ROW, NH, NL ! ! A = REAL ARRAY, PASSED BY COLUMNS ! MAX = MAX NUMBER OF COLUMNS PER SCREEN ! NR = NUMBER OF ROWS IN A ! NC = NUMBER OF COLUMNS IN A ! ! LOOP OVER A SCREEN FULL OF COLUMNS DO J = 1, NC, MAX MAXCOL = 1 K = NC - J - 1 MAXCOL = MIN0 (K,MAX) IF ( NC > MAX ) THEN ! FILL THE TEMPORARY COLUMN NUMBERS TO PRINT DO L = 1, MAXCOL NCOL(L) = L + J - 1 END DO ! OF L COLUMNS PRINT '("ROW/COL", I7, 6I12)', ( NCOL(IC),IC=1,MAXCOL ) END IF ! OF NC VS SCREEN WIDTH ! PRINT ROWS IN CURRENT SCREEN COLUMNS DO ROW = 1, NR NL = ROW + (J-1)*NR NH = NL + (MAXCOL - 1)*NR IF ( NC > MAX ) THEN ! WITH ROW NUMBER WIDER THAN SCREEN PRINT '(I3, 1X, 7(1PG12.5))', ROW,( A(I),I=NL,NH,NR ) ELSE ! NC ! WITHOUT ROW NUMBER FITS SCREEN PRINT '(7(1PG12.5))', ( A(I),I=NL,NH,NR ) END IF ! OF NC VS SCREEN WIDTH END DO ! OF ROW OVER ROWS END DO ! OF J PRINT *, ' ' END SUBROUTINE G_PRINT SUBROUTINE IPRINT (m, NR, NC) ! * * * * * * * * * * * * * * * * * * * * * * * ! PRINTING OF AN INTEGER ARRAY M(NR,NC) ! * * * * * * * * * * * * * * * * * * * * * * * IMPLICIT NONE INTEGER, PARAMETER :: N_PRT = 6, MAX = 10 INTEGER, INTENT(IN) :: NR, NC INTEGER, INTENT(IN) :: m (nr*nc) ! Automatic Arrays INTEGER :: NCOL (MAX) INTEGER :: I, IC, IR, J, L, MAXCOL, NH, NL ! m = INTEGER ARRAY ! NR = NUMBER OF ROWS IN m ! NC = NUMBER OF COLUMNS IN m DO J = 1, NC, MAX MAXCOL = MIN0 ((NC - (J - 1)), MAX) DO L = 1, MAXCOL NCOL (L) = L + (J - 1) END DO WRITE (N_PRT, 5000) (NCOL (IC), IC = 1, MAXCOL) 5000 FORMAT ('ROW/COL',I7, 9I10 ) DO IR = 1, NR NL = IR + (J - 1) * NR NH = NL + (MAXCOL - 1) * NR WRITE (N_PRT, 5010) IR, (m (I), I = NL, NH, NR) 5010 FORMAT ( I4, 10I10 ) END DO END DO END SUBROUTINE IPRINT SUBROUTINE MAX_AND_MIN_RANGE (DD) ! * * * * * * * * * * * * * * * * * * * * * * * * * * ! LIST EXTREME RANGE OF VALUES OF THE N_G_DOF NODAL DOF ! * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants Use Interface_Header IMPLICIT NONE REAL(DP), INTENT(IN) :: DD (N_D_FRE) ! Automatic Arrays REAL(DP) :: RANGE (N_G_DOF, 2), DD_TEST INTEGER :: N_RANGE (N_G_DOF, 2), INDEX (N_G_DOF) INTEGER :: J, I ! MAX_NP = NUMBER OF NODES IN SYSTEM ! N_G_DOF = NUMBER OF PARAMETERS (DOF) PER NODE ! N_D_FRE = TOTAL NUMBER OF SYSTEM DEGREES OF FREEDOM ! RANGE : 1-MAXIMUM VALUE, 2-MINIMUM VALUE ! DD = ARRAY OF SYSTEM DEGREES OF FREEDOM ! INDEX = LIST OF SYSTEM DOF NOS FOR DOF AT NODE ! N_RANGE = ARRAY OF NODE NOS OF EXTREME VALUE POINTS RANGE (1:N_G_DOF, 1) = DD (1:N_G_DOF) ! Initialize RANGE (1:N_G_DOF, 2) = DD (1:N_G_DOF) ! Initialize N_RANGE (1:N_G_DOF, 1) = 1 ! Initialize N_RANGE (1:N_G_DOF, 2) = 1 ! Initialize DO I = 1, MAX_NP INDEX = GET_INDEX_AT_PT (I) DO J = 1, N_G_DOF DD_TEST = DD (INDEX (J) ) IF ( DD_TEST > RANGE (J, 1) ) THEN ! CURRENT MAXIMUM RANGE (J, 1) = DD_TEST ; N_RANGE (J, 1) = I CYCLE ! TO NEXT J ELSE IF ( DD_TEST < RANGE (J, 2) ) THEN ! CURRENT MINIMUM RANGE (J, 2) = DD_TEST ; N_RANGE (J, 2) = I END IF END DO END DO ! PRINT RANGE OF VALUES !b WRITE (U_HIST, 5000) !b WRITE (U_LOG , 5000) WRITE (N_PRT, 5000) ; 5000 FORMAT ( /, & '*** EXTREME VALUES OF THE NODAL PARAMETERS ***',/, & 'PARAMETER MAXIMUM, NODE MINIMUM, NODE') DO J = 1, N_G_DOF !b WRITE (U_HIST, 10) DOF_NAME(J), RANGE (J, 1), N_RANGE (J, 1),& !b RANGE (J, 2), N_RANGE (J, 2) !b WRITE (U_LOG , 10) DOF_NAME(J), RANGE (J, 1), N_RANGE (J, 1),& !b RANGE (J, 2), N_RANGE (J, 2) WRITE (N_PRT, 10) DOF_NAME(J), RANGE (J, 1), N_RANGE (J, 1),& RANGE (J, 2), N_RANGE (J, 2) 10 FORMAT (A12, 2X, 1PE11.4, ',', I7, 2X, 1PE11.4, ',', I7) END DO IF ( .NOT. IS_NULL ) THEN ! Expect non-zero solutions DD_TEST = 2 * TINY (DD_TEST) IF ( MAXVAL (ABS (RANGE)) < DD_TEST ) THEN ! Abort N_WARN = N_WARN + 1 PRINT *, 'WARNING: MAX_AND_MIN_RANGE NULL SOLUTION, MUST' PRINT *, 'DELETE el_no_col, OR ADD is_null IN KEYWORDS' PRINT *, 'ERROR, NULL SOLUTION, CURRENT CONTROLS ARE:' CALL LIST_CONTROLS STOP 'ERROR, NULL SOLUTION. ADD is_null IN KEYWORDS' END IF ! Abort END IF ! not expecting null solution END SUBROUTINE MAX_AND_MIN_RANGE SUBROUTINE MAX_AND_MIN_SCP_AVE_F90 (SCP_AVERAGES) ! * * * * * * * * * * * * * * * * * * * * * * * * * * ! LIST EXTREME RANGE OF VALUES OF THE SCP AVERAGES ! * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants Use Interface_Header IMPLICIT NONE REAL(DP), INTENT(IN) :: SCP_AVERAGES (MAX_NP, SCP_FIT) ! Automatic Arrays REAL(DP) :: RANGE (SCP_FIT, 2), AVE_TEST INTEGER :: N_RANGE (SCP_FIT, 2) INTEGER :: J, I ! MAX_NP = NUMBER OF NODES IN SYSTEM ! SCP_FIT = NUMBER OF AVERAGED FLUX ITEMS PER NODE ! N_D_FRE = TOTAL NUMBER OF SYSTEM DEGREES OF FREEDOM ! RANGE : 1-MAXIMUM VALUE, 2-MINIMUM VALUE ! SCP_AVERAGES = ARRAY OF SYSTEM DEGREES OF FREEDOM ! N_RANGE = ARRAY OF NODE NUMBERS OF EXTREME VALUE POINTS RANGE (1:SCP_FIT, 1) = SCP_AVERAGES (1, 1:SCP_FIT) ! Initialize RANGE (1:SCP_FIT, 2) = SCP_AVERAGES (1, 1:SCP_FIT) ! Initialize N_RANGE (1:SCP_FIT, 1) = 1 ! Initialize N_RANGE (1:SCP_FIT, 2) = 1 ! Initialize DO J = 1, SCP_FIT DO I = 1, MAX_NP AVE_TEST = SCP_AVERAGES (I, J) IF ( AVE_TEST > RANGE (J, 1) ) THEN ! CURRENT MAXIMUM RANGE (J, 1) = AVE_TEST ; N_RANGE (J, 1) = I ELSE IF ( AVE_TEST < RANGE (J, 2) ) THEN ! CURRENT MINIMUM RANGE (J, 2) = AVE_TEST ; N_RANGE (J, 2) = I END IF END DO ! on nodes END DO ! on items ! PRINT RANGE OF VALUES !b WRITE (U_HIST, 5000) !b WRITE (U_LOG , 5000) WRITE (N_PRT, 5000) ; 5000 FORMAT ( /, & '*** EXTREME VALUES OF THE NODAL SCP_AVERAGE ITEMS ***',/, & 'ITEM MAXIMUM, NODE MINIMUM, NODE') DO J = 1, SCP_FIT !b WRITE (U_HIST, 10) AVE_NAME (J), RANGE (J, 1), N_RANGE (J, 1),& !b RANGE (J, 2), N_RANGE (J, 2) !b WRITE (U_LOG , 10) AVE_NAME (J), RANGE (J, 1), N_RANGE (J, 1),& !b RANGE (J, 2), N_RANGE (J, 2) WRITE (N_PRT, 10) AVE_NAME (J), RANGE (J, 1), N_RANGE (J, 1),& RANGE (J, 2), N_RANGE (J, 2) 10 FORMAT (A12, 1PE11.4, ',', I7, 2X, 1PE11.4, ',', I7) END DO AVE_TEST = 2 * TINY (AVE_TEST) ! Small value IF ( MAXVAL (ABS (RANGE)) < AVE_TEST ) THEN ! Warn PRINT *, 'WARNING: MAX_AND_MIN_SCP_AVE, ALL AVERAGES ARE ZERO' PRINT *, ' IN F90 VERSION' PRINT *, ' ERROR ESTIMATOR WILL BE SKIPPED' NULL_SCP_AVE = .TRUE. !b print *, 'NULL_SCP_AVE = ', NULL_SCP_AVE !b SKIP_ERROR = .TRUE. N_WARN = N_WARN + 1 END IF ! Warn END SUBROUTINE MAX_AND_MIN_SCP_AVE_F90 SUBROUTINE MAX_AND_MIN_SCP_AVE_F95 (SCP_AVERAGES) ! * * * * * * * * * * * * * * * * * * * * * * * * * * ! LIST EXTREME RANGE OF VALUES OF THE SCP AVERAGES ! USING F95 EXTENSION OF MAXLOC INTRINSIC FUNCTION ! * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants Use Interface_Header IMPLICIT NONE REAL(DP), INTENT(IN) :: SCP_AVERAGES (MAX_NP, SCP_FIT) REAL(DP) :: MAX_AVE_VALUE (SCP_FIT), MIN_AVE_VALUE (SCP_FIT) INTEGER :: MAX_AVE_NODE (SCP_FIT), MIN_AVE_NODE (SCP_FIT) REAL(DP) :: AVE_TEST INTEGER :: J ! MAX_NP = NUMBER OF NODES IN SYSTEM ! SCP_FIT = NUMBER OF AVERAGED FLUX ITEMS PER NODE ! SCP_AVERAGES = ARRAY OF SYSTEM DEGREES OF FREEDOM MAX_AVE_VALUE = MAXVAL (SCP_AVERAGES, DIM = 1) MAX_AVE_NODE = MAXLOC (SCP_AVERAGES, DIM = 1) MIN_AVE_VALUE = MINVAL (SCP_AVERAGES, DIM = 1) MIN_AVE_NODE = MINLOC (SCP_AVERAGES, DIM = 1) ! PRINT RANGE OF VALUES !b WRITE (U_HIST, 5000) !b WRITE (U_LOG , 5000) WRITE (N_PRT, 5000) ; 5000 FORMAT ( /, & '*** EXTREME VALUES OF THE NODAL SCP_AVERAGE ITEMS ***',/, & 'ITEM MAXIMUM, NODE MINIMUM, NODE') DO J = 1, SCP_FIT !b WRITE (U_HIST, 6) AVE_NAME (J), & !b MAX_AVE_VALUE (J), MAX_AVE_NODE (J), & !b MIN_AVE_VALUE (J), MIN_AVE_NODE (J) !b WRITE (U_LOG , 6) AVE_NAME (J), & !b MAX_AVE_VALUE (J), MAX_AVE_NODE (J), & !b MIN_AVE_VALUE (J), MIN_AVE_NODE (J) WRITE (N_PRT, 6) AVE_NAME (J), & MAX_AVE_VALUE (J), MAX_AVE_NODE (J), & MIN_AVE_VALUE (J), MIN_AVE_NODE (J) 6 FORMAT (A12, 1PE11.4, ',', I7, 2X, 1PE11.4, ',', I7) END DO AVE_TEST = 2 * TINY (AVE_TEST) ! Small value IF ( (ABS (MAXVAL (MAX_AVE_VALUE)) < AVE_TEST) .AND. & (ABS (MINVAL (MIN_AVE_VALUE)) < AVE_TEST) ) THEN ! Warn PRINT *, 'WARNING: MAX_AND_MIN_SCP_AVE, ALL AVERAGES ARE ZERO' PRINT *, ' IN F95 VERSION' PRINT *, ' ERROR ESTIMATOR WILL BE SKIPPED' NULL_SCP_AVE = .TRUE. SKIP_ERROR = .TRUE. N_WARN = N_WARN + 1 END IF ! Warn END SUBROUTINE MAX_AND_MIN_SCP_AVE_F95 SUBROUTINE RPRINT (A, NR, NC, I_OPT) ! * * * * * * * * * * * * * * * * * * * * * * * ! PRINTING OF REAL MATRIX A(NR,NC) ! * * * * * * * * * * * * * * * * * * * * * * * Use Precision_Module IMPLICIT NONE INTEGER, PARAMETER :: N_PRT = 6, MAX = 10 INTEGER, INTENT(IN) :: NR, NC, I_OPT REAL(DP), INTENT(IN) :: A (nr*nc) ! Automatic Arrays INTEGER :: NCOL (MAX) INTEGER :: I, IC, IR, J, L, MAXCOL, NH, NL ! A = REAL ARRAY ! NR = NUMBER OF ROWS IN A ! NC = NUMBER OF COLUMNS IN A !b! I_OPT = 0 USE F FORMAT, OTHERWISE USE E FORMAT ! I_OPT = 0 F FORMAT, 1 E FORMAT, 2 F LONG, 3 E LONG DO J = 1, NC, MAX MAXCOL = MIN0 ((NC - (J - 1)), MAX) DO L = 1, MAXCOL NCOL (L) = L + (J - 1) END DO IF ( I_OPT < 2 ) THEN !b WRITE (N_PRT, 5000) (NCOL (IC), IC = 1, MAXCOL) 5000 FORMAT ('ROW/COL', I7, 9I10 ) ELSE WRITE (N_PRT, 5001) (NCOL (IC), IC = 1, MAXCOL) 5001 FORMAT ('ROW/COL', I11, 9I14 ) END IF DO IR = 1, NR NL = IR + (J - 1) * NR NH = NL + (MAXCOL - 1)* NR IF ( I_OPT == 0 ) THEN WRITE (N_PRT, 5010) IR, (A (I), I = NL, NH, NR) 5010 FORMAT (I4, 10F10.4) !b ELSE ELSEIF ( I_OPT == 1 ) THEN !b WRITE (N_PRT, 5020) IR, (A (I), I = NL, NH, NR) 5020 FORMAT (I4, 10(1PE10.2) ) ELSEIF ( I_OPT == 1 ) THEN !b WRITE (N_PRT, 5011) IR, (A (I), I = NL, NH, NR) !b 5011 FORMAT (I4, 10F14.8) !b ELSEIF ( I_OPT == 3 ) THEN !b WRITE (N_PRT, 5021) IR, (A (I), I = NL, NH, NR) !b 5021 FORMAT (I4, 10(1PE14.6) ) !b END IF END DO END DO END SUBROUTINE RPRINT SUBROUTINE WRITE_BY_ELEMENTS (DD, NODES) ! * * * * * * * * * * * * * * * * * * * * * * * * * * ! OUTPUT, BY ELEMS, OF CALCULATED DEGREES OF FREEDOM ! * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants Use Elem_Type_Data Use Interface_Header IMPLICIT NONE REAL(DP), INTENT(IN) :: DD (N_D_FRE) INTEGER, INTENT(IN) :: NODES (L_S_TOT, NOD_PER_EL) ! Automatic Arrays INTEGER :: INDX (N_G_DOF), EL_NODES (NOD_PER_EL) INTEGER :: IE, K, L, NODE ! N_ELEMS = NUMBER OF ELEMENTS IN SYSTEM ! NOD_PER_EL = NUMBER NODES PER ELEMENT ! N_G_DOF = NUMBER OF PARAMETERS (DOF) PER NODE ! N_D_FRE = NUMBER DEGREES OF FREEDOM IN SYSTEM ! DD = CALCULATED NODAL PARAMETERS (DOF) ! INDX = SYSTEM DOF NUMBERS FOR ELEMENT PARAMETERS ! NODES = NODAL INCIDENCES OF ALL ELEMENTS ! EL_NODES = NODAL INCIDENCES OF AN ELEMENT !b WRITE (N_PRT, 5) N_G_DOF !b 5 FORMAT (/, '*** OUTPUT OF RESULTS BY ELEMENT ORDER ***',/, & !b 'ELEMENT, NODE, ', I3,' PARAMETERS' ) IF ( TRANSIENT ) PRINT *,'At TIME = ', TIME WRITE (N_PRT, 5) DOF_NAME (1:N_G_DOF) 5 FORMAT (/, '*** OUTPUT OF RESULTS BY ELEMENT ORDER ***',/, & 'ELEMENT, NODE, ', (6A12) ) DO IE = 1, N_ELEMS LT_N = LT_DATA (1, L_TYPE(IE)) !9 use function to get EL_NODES = GET_ELEM_NODES (IE, LT_N, NODES) DO K = 1, LT_N NODE = EL_NODES (K) ! ALLOW FOR OMITTED NODES IF ( NODE > 0 ) THEN INDX = GET_INDEX_AT_PT (NODE) WRITE (N_PRT, 6) IE, NODE, (DD (INDX (L) ), L = 1, N_G_DOF) !b 6 FORMAT (I5, I8, 2X, (6(1X, 1PE12.5)) ) 6 FORMAT (I5, I8, 2X, (6(1X, ES11.4)) ) !b END IF END DO END DO END SUBROUTINE WRITE_BY_ELEMENTS SUBROUTINE WRITE_BY_NODES (X, DD) ! * * * * * * * * * * * * * * * * * * * * * * * * * * ! OUTPUT, BY NODES, OF CALCULATED DEGREES OF FREEDOM ! * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants Use Interface_Header IMPLICIT NONE REAL(DP), INTENT(IN) :: X (MAX_NP, N_SPACE), DD (N_D_FRE) INTEGER :: I, K, L ! Automatic Arrays INTEGER :: INDEX (N_G_DOF) ! MAX_NP = NUMBER OF NODES IN SYSTEM ! N_G_DOF = NUMBER OF PARAMETERS (DOF) PER NODE ! N_D_FRE = NUMBER OF DOF IN THE SYSTEM ! N_SPACE = DIMENSION OF SPACE ! X = SYSTEM COORDINATES OF ALL NODES ! DD = CALCULATED NODAL PARAMETERS ! INDEX = SYSTEM DOF NOS OF PARAMETERS ON A NODE IF ( TRANSIENT ) PRINT *,'At TIME = ', TIME IF ( N_SPACE==3 .AND. N_G_DOF == 6 ) THEN ! space frame or shell WRITE (N_PRT, 30) XYZ_NAME(1:3), DOF_NAME(1:3) WRITE (N_PRT, 31) DOF_NAME(4:6) 30 FORMAT( /, '*** OUTPUT OF RESULTS IN NODAL ORDER ***', & /, ' NODE, ', (6A12) ) 31 FORMAT( ' NODE, ', 36X, (3A12) ) DO I = 1, MAX_NP INDEX = GET_INDEX_AT_PT (I) WRITE (N_PRT, 41) I, (X (I, L), L = 1, 3), & (DD (INDEX (K)), K = 1, 3) WRITE (N_PRT, 42) I, (DD (INDEX (K)), K = 4, 6) 41 FORMAT ( I5, (9(1X, ES11.4)) ) !b 42 FORMAT ( I5, 36X, (3(1X, ES11.4)) ) !b END DO ELSE WRITE (N_PRT, 50) XYZ_NAME(1:N_SPACE), DOF_NAME(1:N_G_DOF) 50 FORMAT( /, '*** OUTPUT OF RESULTS IN NODAL ORDER ***', & /, ' NODE, ', (9A12) ) DO I = 1, MAX_NP INDEX = GET_INDEX_AT_PT (I) WRITE (N_PRT, 51) I, (X (I, L), L = 1, N_SPACE), & (DD (INDEX (K)), K = 1, N_G_DOF) 51 FORMAT ( I8, (9(1X, ES11.4)) ) !b END DO END IF ! space frame or shell END SUBROUTINE WRITE_BY_NODES SUBROUTINE WRITE_BY_NODES_WITH_EXACT (X, DD) ! * * * * * * * * * * * * * * * * * * * * * * * * * * ! OUTPUT, BY NODES, OF CALCULATED DEGREES OF FREEDOM ! AND THE EXACT ANSWERS ! * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants Use Select_Source Use Interface_Header IMPLICIT NONE REAL(DP), INTENT(IN) :: X (MAX_NP, N_SPACE), DD (N_D_FRE) INTEGER :: I, K, L ! Automatic Arrays INTEGER :: INDEX (N_G_DOF) REAL(DP) :: EXACT_SOL (N_G_DOF) ! PT EXACT VALUE ! MAX_NP = NUMBER OF NODES IN SYSTEM ! N_G_DOF = NUMBER OF PARAMETERS (DOF) PER NODE ! N_D_FRE = NUMBER OF DOF IN THE SYSTEM ! N_SPACE = DIMENSION OF SPACE ! X = SYSTEM COORDINATES OF ALL NODES ! DD = CALCULATED NODAL PARAMETERS ! INDEX = SYSTEM DOF NOS OF PARAMETERS ON A NODE WRITE (N_PRT, 50) XYZ_NAME (1:N_SPACE), DOF_NAME (1:N_G_DOF), & EXACT_NAME (1:N_G_DOF) ; 50 FORMAT( /, & '*** OUTPUT OF RESULTS AND EXACT VALUES IN NODAL ORDER ***', /, & ' NODE, ', (9A12) ) DO I = 1, MAX_NP CALL SELECT_EXACT_SOLUTION (X (I, :), EXACT_SOL) INDEX = GET_INDEX_AT_PT (I) IF ( N_SPACE < 3 ) THEN WRITE (N_PRT, 51) I, (X (I, L), L = 1, N_SPACE), & (DD (INDEX (K)), K = 1, N_G_DOF), EXACT_SOL ELSE WRITE (N_PRT, 51) I, (X (I, L), L = 1, N_SPACE), & (DD (INDEX (K)), K = 1, N_G_DOF) WRITE (N_PRT, 52) I, EXACT_SOL END IF 51 FORMAT ( I6, (9(1X, ES11.4)) ) !b 52 FORMAT ( I6, 36X, (6(1X, ES11.4)) ) !b END DO END SUBROUTINE WRITE_BY_NODES_WITH_EXACT SUBROUTINE DISPLAY_D_MATRIX ( A ) ! ------------------------------------------------------ ! PRETTY PRINTING OF DP MATRIX A(NR,NC) WITH E FORMAT ! ------------------------------------------------------ Use Precision_Module IMPLICIT NONE REAL(DP), INTENT(IN) :: A(:,:) ! THE ARRAY INTEGER, PARAMETER :: MAX_COL = 7 ! COLS PER PAGE INTEGER :: NCOL(MAX_COL) ! COLUMN NUMBERS INTEGER :: NC, NR, COL, N, J, K, FULL INTEGER :: PART, COL_SIZE, PAGES ! A = REAL ARRAY, PASSED BY COLUMNS ! COL = MAX NUMBER OF COLUMNS PER SCREEN ! NC = NUMBER OF COLUMNS IN A ! NCOL = VECTOR SUBSCRIPT FOR COLUMNS TO BE PRINTED ! NR = NUMBER OF ROWS IN A NR = SIZE (A, 1) ! FIND NUMBER OF ROWS NC = SIZE (A, 2) ! FIND NUMBER OF COLUMNS nr = SIZE (A, DIM=1) NC = SIZE (A, DIM=2) if ( NR < 1 .OR. NC < 1 ) STOP 'DISPLAY_D_MATRIX: zero size' COL = MIN0(NC,MAX_COL) ! LOOP OVER A SCREEN FULL OF COLUMNS PAGES = CEILING(REAL(NC)/COL) FULL = NC/COL PART = NC - FULL*COL DO J = 1, PAGES ! FILL THE TEMPORARY COLUMN NUMBERS TO PRINT IF ( J > FULL ) THEN ! LAST PARTIAL PAGE COL_SIZE = PART ELSE ! FULL PAGE COL_SIZE = COL END IF ! STATUS CHECK !b NCOL = (/ (((J-1)*COL+K), K = 1,COL_SIZE) /) DO K = 1, COL_SIZE NCOL (K) = (J-1)*COL + K END DO WRITE (*,5) ( NCOL(K), K = 1,COL_SIZE) 5 FORMAT (' R\C', I6, 11I12 ) ! PRINT ROWS IN CURRENT SCREEN COLUMNS DO N = 1, NR WRITE (*,10) N, ( A(N,NCOL(K)), K = 1,COL_SIZE ) !b 10 FORMAT ( I3, 11(1PE12.4) ) 10 FORMAT ( I3, 11(ES12.4) ) !b END DO ! OF N OVER ROWS END DO ! OF J COLUMNS WRITE (*,*) ' ' END SUBROUTINE DISPLAY_D_MATRIX SUBROUTINE DISPLAY_D_VECTOR ( A ) ! ------------------------------------------------------ ! PRETTY PRINTING OF DP VECTOR A(NC) WITH E FORMAT ! ------------------------------------------------------ Use Precision_Module IMPLICIT NONE REAL(DP), INTENT(IN) :: A(:) ! THE ARRAY INTEGER, PARAMETER :: MAX_COL = 7 ! COLS PER PAGE INTEGER :: NCOL(MAX_COL) ! COLUMN NUMBERS INTEGER :: NC, MAX, N=1, J, K, FULL INTEGER :: PART, COL_SIZE, PAGES ! A = REAL ARRAY, PASSED BY COLUMNS ! MAX = MAX NUMBER OF COLUMNS PER SCREEN ! NC = NUMBER OF COLUMNS IN A ! NCOL = VECTOR SUBSCRIPT FOR COLUMNS TO BE PRINTED NC = SIZE (A) ! FIND NUMBER OF ROWS if ( NC < 1 ) STOP 'DISPLAY_D_VECTOR: zero size' MAX = MIN0(NC,MAX_COL) ! LOOP OVER A SCREEN FULL OF COLUMNS PAGES = CEILING(REAL(NC)/MAX) FULL = NC/MAX PART = NC - FULL*MAX DO J = 1, PAGES ! FILL THE TEMPORARY COLUMN NUMBERS TO PRINT IF ( J > FULL ) THEN ! LAST PARTIAL PAGE COL_SIZE = PART ELSE ! FULL PAGE COL_SIZE = MAX END IF ! STATUS CHECK !b NCOL = (/ ((J-1)*MAX+K, K = 1,COL_SIZE) /) DO K = 1, COL_SIZE NCOL (K) = (J-1)*MAX + K END DO WRITE (*,5) ( NCOL(K), K = 1,COL_SIZE) 5 FORMAT (' R\C', I6, 11I12 ) ! PRINT ROW IN CURRENT SCREEN COLUMNS WRITE (*,10) N, ( A(NCOL(K)), K = 1,COL_SIZE ) 10 FORMAT ( I3, 11(1PE12.4) ) END DO ! OF J COLUMNS WRITE (*,*) ' ' END SUBROUTINE DISPLAY_D_VECTOR SUBROUTINE DISPLAY_R_MATRIX ( A ) ! ------------------------------------------------------ ! PRETTY PRINTING OF REAL MATRIX A(NR,NC) WITH E FORMAT ! ------------------------------------------------------ IMPLICIT NONE REAL, INTENT(IN) :: A(:,:) ! THE ARRAY INTEGER, PARAMETER :: MAX_COL = 7 ! COLS PER PAGE INTEGER :: NCOL(MAX_COL) ! COLUMN NUMBERS INTEGER :: NC, NR, MAX, N, J, K, FULL INTEGER :: PART, COL_SIZE, PAGES ! A = REAL ARRAY, PASSED BY COLUMNS ! MAX = MAX NUMBER OF COLUMNS PER SCREEN ! NC = NUMBER OF COLUMNS IN A ! NCOL = VECTOR SUBSCRIPT FOR COLUMNS TO BE PRINTED ! NR = NUMBER OF ROWS IN A NR = SIZE (A, 1) ! FIND NUMBER OF ROWS NC = SIZE (A, 2) ! FIND NUMBER OF COLUMNS MAX = MIN0(NC,MAX_COL) ! LOOP OVER A SCREEN FULL OF COLUMNS PAGES = CEILING(REAL(NC)/MAX) FULL = NC/MAX PART = NC - FULL*MAX DO J = 1, PAGES ! FILL THE TEMPORARY COLUMN NUMBERS TO PRINT IF ( J > FULL ) THEN ! LAST PARTIAL PAGE COL_SIZE = PART ELSE ! FULL PAGE COL_SIZE = MAX END IF ! STATUS CHECK !b NCOL = (/ ((J-1)*MAX+K, K = 1,COL_SIZE) /) DO K = 1, COL_SIZE NCOL (K) = (J-1)*MAX + K END DO WRITE (*,5) ( NCOL(K), K = 1,COL_SIZE) 5 FORMAT (' R\C', I6, 11I12 ) ! PRINT ROWS IN CURRENT SCREEN COLUMNS DO N = 1, NR WRITE (*,10) N, ( A(N,NCOL(K)), K = 1,COL_SIZE ) 10 FORMAT ( I3, 11(1PE12.4) ) END DO ! OF N OVER ROWS END DO ! OF J COLUMNS WRITE (*,*) ' ' END SUBROUTINE DISPLAY_R_MATRIX SUBROUTINE DISPLAY_R_VECTOR ( A ) ! ------------------------------------------------------ ! PRETTY PRINTING OF REAL MATRIX A(NR,NC) WITH E FORMAT ! ------------------------------------------------------ IMPLICIT NONE REAL, INTENT(IN) :: A(:) ! THE ARRAY INTEGER, PARAMETER :: MAX_COL = 7 ! COLS PER PAGE INTEGER :: NCOL(MAX_COL) ! COLUMN NUMBERS INTEGER :: NC, MAX, N=1, J, K, FULL INTEGER :: PART, COL_SIZE, PAGES ! A = REAL ARRAY, PASSED BY COLUMNS ! MAX = MAX NUMBER OF COLUMNS PER SCREEN ! NC = NUMBER OF COLUMNS IN A ! NCOL = VECTOR SUBSCRIPT FOR COLUMNS TO BE PRINTED NC = SIZE (A) ! FIND NUMBER OF ROWS MAX = MIN0(NC,MAX_COL) ! LOOP OVER A SCREEN FULL OF COLUMNS PAGES = CEILING(REAL(NC)/MAX) FULL = NC/MAX PART = NC - FULL*MAX DO J = 1, PAGES ! FILL THE TEMPORARY COLUMN NUMBERS TO PRINT IF ( J > FULL ) THEN ! LAST PARTIAL PAGE COL_SIZE = PART ELSE ! FULL PAGE COL_SIZE = MAX END IF ! STATUS CHECK !b NCOL = (/ ((J-1)*MAX+K, K = 1,COL_SIZE) /) DO K = 1, COL_SIZE NCOL (K) = (J-1)*MAX + K END DO WRITE (*,5) ( NCOL(K), K = 1,COL_SIZE) 5 FORMAT (' R\C', I6, 11I12 ) ! PRINT ROW IN CURRENT SCREEN COLUMNS WRITE (*,10) N, ( A(NCOL(K)), K = 1,COL_SIZE ) 10 FORMAT ( I3, 11(1PE12.4) ) END DO ! OF J COLUMNS WRITE (*,*) ' ' END SUBROUTINE DISPLAY_R_VECTOR SUBROUTINE DISPLAY_I_MATRIX ( A ) ! ------------------------------------------------------ ! PRETTY PRINTING OF INTEGER MATRIX A(NR,NC) WITH I FORMAT ! ------------------------------------------------------ IMPLICIT NONE INTEGER, INTENT(IN) :: A(:,:) ! THE ARRAY INTEGER, PARAMETER :: MAX_COL = 7 ! COLS PER PAGE INTEGER, ALLOCATABLE :: NCOL(:) ! COLUMN NUMBERS INTEGER NC, NR, MAX, N, J, K, FULL, PART, COL_SIZE, PAGES ! A = INTEGER ARRAY ! MAX = MAX NUMBER OF COLUMNS PER SCREEN ! NC = NUMBER OF COLUMNS IN A ! NCOL = VECTOR SUBSCRIPT FOR COLUMNS TO BE PRINTED ! NR = NUMBER OF ROWS IN A NR = SIZE (A, 1) ! FIND NUMBER OF ROWS NC = SIZE (A, 2) ! FIND NUMBER OF COLUMNS MAX = MIN0(NC,MAX_COL) ALLOCATE( NCOL(MAX) ) ! LOOP OVER A SCREEN FULL OF COLUMNS PAGES = CEILING(REAL(NC)/MAX) FULL = NC/MAX PART = NC - FULL*MAX DO J = 1, PAGES ! FILL THE TEMPORARY COLUMN NUMBERS TO PRINT IF ( J > FULL ) THEN ! LAST PARTIAL PAGE COL_SIZE = PART ELSE ! FULL PAGE COL_SIZE = MAX END IF ! STATUS CHECK !b NCOL = (/ ((J-1)*MAX+K, K = 1,COL_SIZE) /) DO K = 1, COL_SIZE NCOL (K) = (J-1)*MAX + K END DO WRITE (*,5) ( NCOL(K), K = 1,COL_SIZE) 5 FORMAT (' R\C', I6, 11I11 ) ! PRINT ROWS IN CURRENT SCREEN COLUMNS DO N = 1, NR WRITE (*,10) N, ( A(N,NCOL(K)), K = 1,COL_SIZE ) 10 FORMAT ( I3, 11( I11 ) ) END DO ! OF N OVER ROWS END DO ! OF J COLUMNS WRITE (*,*) ' ' END SUBROUTINE DISPLAY_I_MATRIX SUBROUTINE DISPLAY_I_VECTOR ( A ) ! ------------------------------------------------------ ! PRETTY PRINTING OF INTEGER VECTOR A(NC) WITH I FORMAT ! ------------------------------------------------------ IMPLICIT NONE INTEGER, INTENT(IN) :: A(:) ! THE ARRAY INTEGER, PARAMETER :: MAX_COL = 7 ! COLS PER PAGE INTEGER, ALLOCATABLE :: NCOL(:) ! COLUMN NUMBERS INTEGER :: NC, MAX, N=1, J, K, FULL, PART, COL_SIZE, PAGES ! A = INTEGER ARRAY ! MAX = MAX NUMBER OF COLUMNS PER SCREEN ! NC = NUMBER OF COLUMNS IN A ! NCOL = VECTOR SUBSCRIPT FOR COLUMNS TO BE PRINTED NC = SIZE (A) ! FIND NUMBER OF ROWS MAX = MIN0(NC,MAX_COL) ALLOCATE( NCOL(MAX) ) ! LOOP OVER A SCREEN FULL OF COLUMNS PAGES = CEILING(REAL(NC)/MAX) FULL = NC/MAX PART = NC - FULL*MAX DO J = 1, PAGES ! FILL THE TEMPORARY COLUMN NUMBERS TO PRINT IF ( J > FULL ) THEN ! LAST PARTIAL PAGE COL_SIZE = PART ELSE ! FULL PAGE COL_SIZE = MAX END IF ! STATUS CHECK !b NCOL = (/ ((J-1)*MAX+K, K = 1,COL_SIZE) /) DO K = 1, COL_SIZE NCOL (K) = (J-1)*MAX + K END DO WRITE (*,5) ( NCOL(K), K = 1,COL_SIZE) 5 FORMAT (' R\C', I6, 11I11 ) ! PRINT ROWS IN CURRENT SCREEN COLUMNS WRITE (*,10) N, ( A(NCOL(K)), K = 1,COL_SIZE ) 10 FORMAT ( I3, 11( I11 ) ) END DO ! OF J COLUMNS WRITE (*,*) ' ' END SUBROUTINE DISPLAY_I_VECTOR SUBROUTINE ELEM_BAR_CHART (I_BAR, I_PARM, G_PARM, PARM) ! * * * * * * * * * * * * * * * * * * * * * * * * * * ! PRINT-PLOT BAR CHARTS OF ELEM PARAMETER I_PARM ! FROM A TOTAL OF G_PARM PER ELEMENT ! * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants IMPLICIT NONE INTEGER, PARAMETER :: MID = 25, LINE = 2 * MID INTEGER, INTENT(IN) :: I_BAR, I_PARM, G_PARM REAL(DP), INTENT(IN) :: PARM (N_ELEMS * G_PARM) ! Automatic Arrays INTEGER :: ELEM_BAR (I_BAR) ! Locals INTEGER :: I, INDEX, IPMNG, ISPACE, JY, JZ, & K, KOUNT, L, LIMIT, NLAST, NOPLOT REAL(DP) :: DTEST, RANGE, YMAX, YMIN CHARACTER ALINE (LINE+1), SKIP (LINE+1) CHARACTER(1), SAVE :: BLANK = ' ', DOT = '+', DASH = '-', & A_X = 'X', PLUS = '+', A_0 = '0' DATA NOPLOT / 0 / ! I_PARM = NODAL PARAMETER TO BE GRAPHED,1<=I_PARM<=G_PARM ! ELEM_BAR = LIST OF ELEMENTS TO BE USED (WHEN I_BAR > 1) ! PARM = ARRAY OF ALL NODAL PARAMETERS IN THE SYSTEM ! G_PARM = NUMBER OF PARAMETERS PER ELEMENT ! N_SPACE = DIMENSION OF SOLUTION SPACE ! I_BAR = NUMBER OF ELEMENTS TO BE INCLUDED IN BAR CHART ! IF I_BAR=1 USE ALL, ELEM_BAR NOT USED NOPLOT = NOPLOT + 1 LIMIT = I_BAR IF ( I_BAR == 1) LIMIT = N_ELEMS WRITE (N_PRT, 5000) NOPLOT, I_PARM, LIMIT 5000 FORMAT ('*** PRINT PLOT NUMBER',I3,' ***',/, & 'ELEMENT PARAMETER', I3, & ', EVALUATED AT ',I5,' ELEMENTS',/) IPMNG = I_PARM - G_PARM !--> ESTABLISH GRAPH LIMITING VALUES ! INITALIZE MAX MIN VALUES OF PARAMETER L = 1 IF ( I_BAR /= 1) L = ELEM_BAR (1) INDEX = G_PARM * L + IPMNG YMAX = PARM (INDEX) YMIN = YMAX DO I = 1, LIMIT L = I IF ( I_BAR /= 1) L = ELEM_BAR (I) INDEX = G_PARM * L + IPMNG DTEST = PARM (INDEX) IF ( DTEST > YMAX) YMAX = DTEST IF ( DTEST < YMIN) YMIN = DTEST END DO KOUNT = 1 WRITE (N_PRT, 5010) YMIN, YMAX 5010 FORMAT (' RANGE ON GRAPH IS ',1PE12.5,' TO ',1PE12.5/) IF ( YMIN == YMAX ) THEN WRITE (N_PRT, 5011) 5011 FORMAT ('CONSTANT VALUE, PLOT SKIPPED') RETURN END IF RANGE = YMAX - YMIN ! FLAG ZERO VALUE IF ( YMIN < 0.0 .AND. YMAX > 0.0 ) THEN JZ = 1*(YMAX - 0.0) / RANGE & + (LINE+1)*(0.0 - YMIN) / RANGE + 0.1 ELSE JZ = 0 END IF ! SCALING COMPLETE SKIP (2:LINE) = BLANK ALINE (2:LINE) = DASH SKIP (1) = PLUS SKIP (MID+1) = DOT SKIP (LINE+1) = PLUS ALINE (1:LINE+1:5) = PLUS ALINE (MID+1) = DOT WRITE (N_PRT, 5020) ALINE 5020 FORMAT(' ELEMENT VALUE ', (101A1) ) L = 1 IF ( I_BAR > 1) L = ELEM_BAR (1) DO K = 1, LIMIT NLAST = L L = K IF ( I_BAR > 1) L = ELEM_BAR (K) INDEX = G_PARM * L + IPMNG DTEST = PARM (INDEX) ALINE (2:LINE) = BLANK ALINE (1) = PLUS ALINE (LINE+1) = PLUS IF ( JZ > 0) ALINE (JZ) = A_0 ISPACE = 1 ! LINEAR INTERPOLATION FOR COLUMN NUMBER JY = 1 * (YMAX - DTEST) / RANGE & + (LINE+1) * (DTEST - YMIN) / RANGE + 0.1 ALINE (1:JY) = A_X WRITE (N_PRT, 5040) L, PARM (INDEX), ALINE 5040 FORMAT (I7,2X,1PE10.3,1X,101A1) END DO ALINE (2:LINE) = DASH ALINE (1:LINE+1:5) = PLUS ALINE (MID+1) = DOT WRITE (N_PRT, 5020) ALINE END SUBROUTINE ELEM_BAR_CHART ! now in scp_lib SAVE_QP_FLUX_FOR_SCP (IE, IQ, XYZ, E, STRAIN) SUBROUTINE LIST_ELEM_GRADIENTS (N_FILE, IE) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! LIST ELEMENT GRADIENTS AT QUADRATURE POINTS, ON N_FILE ! GENERALLY USED ONLY FOR ELEM POST_PROCESSING SINCE ! FLUXES ARE NOW STANDARD OUTPUTS ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants Use Elem_Type_Data ! for COORD (LT_N, N_SPACE), D (LT_FREE) IMPLICIT NONE ! ALWAYS USED INTEGER, INTENT(IN) :: N_FILE, IE ! Automatic Arrays REAL(DP) :: XYZ (N_SPACE), DGH (N_SPACE, LT_FREE), & STRAIN (N_R_B + 2), E (N_R_B, N_R_B), RMS REAL(DP) :: B (N_R_B, LT_FREE) !b INTEGER, SAVE :: TEST_IE, TEST_IP, J, N_IP, EOF, IO_1 REAL(DP), SAVE :: DERIV_MAX = 0.d0 ! VARIABLES: ! B = GRADIENT VERSUS DOF MATRIX (N_R_B, LT_FREE) ! COORD = SPATIAL COORDINATES OF ELEMENT'S NODES ! D = NODAL PARAMETERS ASSOCIATED WITH AN ELEMENT ! DGH = GLOBAL DERIVATIVES INTERPOLATION FUNCTIONS ! E = ELEMENT CONSTITUTIVE MATRIX AT GAUSS POINT ! G = GEOMETRIC INTERPOLATION FUNCTIONS ! H = SOLUTION INTERPOLATION FUNCTIONS ! LT_N = NUMBER OF NODES PER ELEMENT ! LT_FREE = NUMBER OF DEGREES OF FREEDOM PER ELEMENT ! N_ELEMS = TOTAL NUMBER OF ELEMENTS ! N_R_B = NUMBER OF ROWS IN B AND E MATRICES ! N_SPACE = DIMENSION OF SPACE ! N_FILE = UNIT FOR POST SOLUTION MATRICES STORAGE ! STRAIN = STRAIN OR GRADIENT VECTOR ! XYZ = SPACE COORDINATES AT A POINT ! U_PLT4 = UNIT TO STORE PLOT DATA, IF > 0 !--> PRINT TITLES ON THE FIRST CALL IF ( IE == 1) THEN REWIND (N_FILE) IF ( LIST_QP_FLUX ) THEN WRITE (N_PRT, 5) N_SPACE, N_R_B 5 FORMAT (/,'*** GRADIENTS AT ELEMENT INTEGRATION POINTS ***',/, & ' ELEMENT, POINT, ', I2, ' COORDINATES', I2, ' GRADIENTS', /) END IF ! OPEN PLOT FILE IF ACTIVE IF ( U_PLT4 > 0 ) OPEN (U_PLT4, FILE='el_qp_xyz_grads.tmp', & ACTION='WRITE', STATUS='REPLACE', IOSTAT = IO_1) IF ( IO_1 > 0 ) PRINT *, 'OPEN FAILED FOR el_qp_xyz_grads.tmp' ENDIF ! IS THIS AN ELEMENT, OR A BOUNDARY SEGMENT ? IF ( IE <= N_ELEMS ) THEN ! ELEMENT RESULTS ! READ NUMBER OF POINTS TO BE CALCULATED READ (N_FILE, IOSTAT = EOF) N_IP IF ( EOF /= 0 ) THEN PRINT *, 'LIST_ELEM_GRADIENTS EOF AT ELEMENT ', IE, EOF PRINT *, 'WAS EXPECTING: N_IP. WAS UNIT ', N_FILE, ' REWOUND ?' STOP 'ERROR, EOF IN LIST_ELEM_GRADIENTS' END IF !--> READ COORDS. AND DERIVATIVE MATRIX DO J = 1, N_IP READ (N_FILE, IOSTAT = EOF) XYZ, E, B !b DGH IF ( EOF /= 0 ) THEN PRINT *, 'LIST_ELEM_GRADIENTS EOF AT ELEMENT ', IE, EOF PRINT *, 'WAS EXPECTING: XYZ, E, B AT POINT ', J, ' OF ', N_IP !b PRINT *, 'WAS EXPECTING: XYZ, E, DGH' CALL INQUIRE_ABOUT_UNIT (N_FILE) STOP 'ERROR, EOF IN LIST_ELEM_GRADIENTS' END IF ! CALCULATE DERIVATIVES, STRAIN = B*D STRAIN (1:N_R_B) = MATMUL (B, D) ! SCALAR MEASURE (VALID FOR VECTORS) RMS = SQRT ( SUM (STRAIN (1:N_R_B)**2) ) !--> PRINT COORDINATES AND GRADIENT AT THE POINT IF ( LIST_QP_FLUX ) THEN WRITE (N_PRT, 20) IE, J, XYZ, STRAIN (1:N_R_B) 20 FORMAT ( I5, I3, 10(1PE11.3), (8X, 10(1PE11.3))) END IF IF ( RMS > DERIV_MAX ) THEN DERIV_MAX = RMS TEST_IE = IE ; TEST_IP = J END IF !--> STORE RESULTS TO BE PLOTTED LATER IF ( U_PLT4 > 0 .AND. IO_1 == 0 ) & WRITE (U_PLT4, 40) XYZ, STRAIN (1:N_R_B) 40 FORMAT ( (10(1PE13.5)) ) END DO PRINT *, ' ' !--> ARE GRADIENT CALCULATIONS COMPLETE FOR ALL ELEMENTS IF (IE == N_ELEMS) THEN ! LIST LARGEST GRADIENT WRITE (N_PRT, "('LARGEST GRADIENT = ', 1PE13.5)") DERIV_MAX WRITE (N_PRT, "('AT ELEM =', I6, ', POINT = ', I2)") & TEST_IE, TEST_IP MAX_DERIVATIVE = DERIV_MAX ! plot scale IF ( U_PLT4 > 0 ) CLOSE (U_PLT4) END IF ! LAST ELEMENT END IF ! ELEMENT OR BOUNDARY SEGMENT END SUBROUTINE LIST_ELEM_GRADIENTS SUBROUTINE LIST_ELEM_SOLUTION_INTEGRAL (N_FILE, IE, EACH) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! LIST INTEGRAL OF SOLUTION FROM DATA POINTS, ON N_FILE ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants Use Elem_Type_Data ! for D (LT_FREE) IMPLICIT NONE ! ALWAYS USED INTEGER, INTENT(IN) :: N_FILE, IE LOGICAL, INTENT(IN) :: EACH ! list each if true REAL(DP), SAVE :: VALUE, TOTAL REAL(DP) :: H_INTG (LT_FREE) INTEGER :: EOF ! VARIABLES: ! D = NODAL PARAMETERS ASSOCIATED WITH AN ELEMENT ! H = SOLUTION INTERPOLATION FUNCTIONS ! H_INTG = INTEGRAL OF INTERPOLATION FUNCTIONS ! LT_FREE = NUMBER OF DEGREES OF FREEDOM PER ELEMENT ! N_FILE = UNIT FOR POST SOLUTION MATRICES STORAGE !--> PRINT TITLES ON THE FIRST CALL AND INITIALIZE IF ( IE == 1) THEN REWIND (N_FILE) TOTAL = 0.d0 IF ( EACH ) WRITE (N_PRT, 5) 5 FORMAT (/,'** VOLUME INTEGRAL OF DEPENDENT VARIABLE **',/, & 'ELEMENT INTEGRAL') ENDIF IF ( IE <= N_ELEMS ) THEN ! ELEMENT RESULTS !--> CALCULATE ELEMENT CONTRIBUTION READ (N_FILE, IOSTAT = EOF) H_INTG IF ( EOF /= 0 ) THEN PRINT *, 'LIST_ELEM_SOLUTION_INTEGRAL EOF AT ELEMENT ', IE PRINT *, 'WAS EXPECTING: H_INTG' STOP 'ERROR, EOF IN LIST_ELEM_SOLUTION_INTEGRAL' END IF VALUE = DOT_PRODUCT (H_INTG, D) TOTAL = TOTAL + VALUE IF ( EACH ) WRITE (N_PRT, '(I7, 1PE18.6)') IE, VALUE IF (IE == N_ELEMS) WRITE (N_PRT, & "('TOTAL SOLUTION VOLUME INTEGRAL = ', 1PE16.6, /)") TOTAL END IF END SUBROUTINE LIST_ELEM_SOLUTION_INTEGRAL SUBROUTINE LIST_ELEM_FLUXES (N_FILE, IE) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! LIST ELEMENT FLUXES AT QUADRATURE POINTS, ON N_FILE ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants Use Elem_Type_Data ! for LT_FREE, D (LT_FREE) IMPLICIT NONE ! ALWAYS USED INTEGER, INTENT(IN) :: N_FILE, IE ! Automatic Arrays !b REAL(DP) :: DGH (N_SPACE, LT_FREE), STRAIN (N_R_B + 2), & REAL(DP) :: STRAIN (N_R_B + 2), & XYZ (N_SPACE), E (N_R_B, N_R_B), STRESS (N_R_B + 2) REAL(DP) :: B (N_R_B, LT_FREE), RMS !b INTEGER, SAVE :: TEST_IE, TEST_IP, J, N_IP, EOF, IO_1 REAL(DP), SAVE :: DERIV_MAX = 0.d0 ! VARIABLES: ! D = NODAL PARAMETERS ASSOCIATED WITH AN ELEMENT ! E = ELEMENT CONSTITUTIVE MATRIX AT GAUSS POINT ! G = GEOMETRIC INTERPOLATION FUNCTIONS ! H = SOLUTION INTERPOLATION FUNCTIONS ! LT_N = NUMBER OF NODES PER ELEMENT ! LT_FREE = NUMBER OF DEGREES OF FREEDOM PER ELEMENT ! N_ELEMS = TOTAL NUMBER OF ELEMENTS ! N_R_B = NUMBER OF ROWS IN B AND E MATRICES ! N_SPACE = DIMENSION OF SPACE ! N_FILE = UNIT FOR POST SOLUTION MATRICES STORAGE ! RMS = ROOT MEAN SQUARE VALUE ! STRAIN = STRAIN OR FLUX VECTOR ! XYZ = SPACE COORDINATES AT A POINT ! U_PLT4 = UNIT TO STORE PLOT DATA, IF > 0 !--> PRINT TITLES ON THE FIRST CALL IF ( IE == 1) THEN REWIND (N_FILE) IF ( LIST_QP_FLUX ) THEN WRITE (N_PRT, 5) XYZ_NAME (1:N_SPACE), FLUX_NAME (1:N_R_B) !b 5 FORMAT (/, & '*** FLUX COMPONENTS AT ELEMENT INTEGRATION POINTS ***',/, & '(SAVED FOR AVERAGING OR APPLICATION POST PROCESSING)',/, & 'ELEMENT, PT, ', (6A12) ) END IF RECORD_NUMBER = 0 ! OPEN FLUX PLOT FILE IF ACTIVE IF ( N_FILE5 > 0 ) OPEN (N_FILE5, FILE='el_qp_xyz_grads.tmp', & ACTION='WRITE', STATUS='REPLACE', IOSTAT = IO_1) IF ( U_PLT4 > 0 ) OPEN (U_PLT4, FILE='el_qp_xyz_fluxes.tmp', & ACTION='WRITE', STATUS='REPLACE', IOSTAT = IO_1) IF ( IO_1 > 0 ) PRINT *, 'OPEN FAILED FOR el_qp_xyz_fluxes.tmp' ENDIF !b print *,'N_FILE, IO_1 ', N_FILE, IO_1 !b !b call LIST_SYSTEM_UNIT_NUMBERS !b ! IS THIS AN ELEMENT, OR A BOUNDARY SEGMENT ? IF ( IE <= N_ELEMS ) THEN ! ELEMENT RESULTS, or IF (IS_ELEMENT) ! READ NUMBER OF INTEGRATION POINTS TO BE EMPLOYED READ (N_FILE, IOSTAT = EOF) N_IP ! expect N_FILE == U_FLUX IF ( EOF /= 0 ) THEN ! HIT END OF FILE PRINT *, 'LIST_ELEM_FLUXES EOF AT ELEMENT ', IE, EOF PRINT *, 'WAS EXPECTING: N_IP. WAS UNIT ', N_FILE, ' REWOUND ?' CALL INQUIRE_ABOUT_UNIT (N_FILE) STOP 'ERROR, EOF IN LIST_ELEM_FLUXES' END IF IF ( DEBUG_SCP ) THEN PRINT *,'ELEM, DOF ', IE, D END IF !--> READ COORDS, CONSTITUTIVE, AND DERIVATIVE MATRIX DO J = 1, N_IP ! OVER ALL INTEGRATION POINTS READ (N_FILE, IOSTAT = EOF) XYZ, E, B ! expect N_FILE == U_FLUX IF ( EOF /= 0 ) THEN PRINT *, 'LIST_ELEM_FLUXES EOF AT ELEMENT ', IE, EOF PRINT *, 'WAS EXPECTING: XYZ, E, B AT POINT ', J, ' OF ', N_IP PRINT *, 'AT UNIT ', N_FILE CALL INQUIRE_ABOUT_UNIT (N_FILE) STOP 'ERROR, EOF IN LIST_ELEM_FLUXES' END IF ! CALCULATE DERIVATIVES, STRAIN = B * D STRAIN (1:N_R_B) = MATMUL (B, D) !b IF ( DEBUG_SCP ) THEN PRINT *,'J, XYZ ', J, XYZ PRINT *,'B MATRIX'; CALL RPRINT (B, N_R_B, LT_FREE, 1) PRINT *,'E MATRIX'; CALL RPRINT (E, N_R_B, N_R_B, 1) PRINT *,'STRAINS '; CALL RPRINT (STRAIN(1:N_R_B), 1, N_R_B, 1) END IF ! FLUX FROM CONSTITUTIVE DATA STRESS (1:N_R_B) = MATMUL (E, STRAIN (1:N_R_B)) !bb IF ( GRAD_BASE_ERROR ) STRESS (1:N_R_B) = STRAIN (1:N_R_B) !b IF ( GRAD_BASE_ERROR ) THEN !b STRESS (1:N_R_B) = STRAIN (1:N_R_B) !b ELSE !b STRESS (1:N_R_B) = MATMUL (E, STRAIN (1:N_R_B)) !b END IF ! GRAD VS FLUX ! SCALAR MEASURE (VALID FOR VECTORS) RMS = SQRT ( SUM (STRESS (1:N_R_B)**2) ) !--> PRINT COORDINATES AND FLUX AT THE POINT IF ( LIST_QP_FLUX ) WRITE (N_PRT, '(I7, I3, 10(ES12.4))') & IE, J, XYZ, STRESS (1:N_R_B) IF ( RMS > DERIV_MAX ) THEN DERIV_MAX = RMS ! LARGEST VALUE TEST_IE = IE ; TEST_IP = J ! AT ELEM AND POINT END IF !--> STORE FLUX RESULTS TO BE PLOTTED LATER, IF USED !b print *,'U_PLT4, IO_1 ', U_PLT4, IO_1 !b IF ( U_PLT4 > 0 .AND. IO_1 == 0 ) & WRITE (U_PLT4, '( (10(1PE13.5)) )') XYZ, STRESS (1:N_R_B) IF ( N_FILE5 > 0 .AND. IO_1 == 0 ) & WRITE (N_FILE5, '( (10(1PE13.5)) )') XYZ, STRAIN (1:N_R_B) ! SAVE COORDINATES & FLUX COMPONENTS FOR SCP FLUX AVERAGING IF ( U_SCPR > 0 ) THEN RECORD_NUMBER = RECORD_NUMBER + 1 SCP_RECORD_NUMBER (IE, J) = RECORD_NUMBER !b WRITE (U_SCPR, REC = RECORD_NUMBER) XYZ, STRESS (1:SCP_FIT) IF ( GRAD_BASE_ERROR ) THEN !b JEA 4/6/02 XXX WRITE (U_SCPR, REC = RECORD_NUMBER) XYZ, STRAIN (1:SCP_FIT) ELSE WRITE (U_SCPR, REC = RECORD_NUMBER) XYZ, STRESS (1:SCP_FIT) END IF ! GRAD VS FLUX ! NOTE: THE ERROR ESTIMATOR WILL NEED E ALSO !b WRITE (U_SCPR, REC = RECORD_NUMBER) XYZ, E, STRESS (1:SCP_FIT) !b change SCP to read E, and change record length !b END IF ! SCP RECOVERY END DO ! OVER INTEGRATION POINTS !b print *,'rec numb ', RECORD_NUMBER !b PRINT *, ' ' !b IF ( I_BUG > 0 ) CALL INQUIRE_ABOUT_SCP_UNIT !b !--> ARE FLUX CALCULATIONS COMPLETE FOR ALL ELEMENTS IF ( IE == N_ELEMS ) THEN ! LIST LARGEST FLUX WRITE (N_PRT, "('LARGEST FLUX RMS_VALUE = ', 1PE13.5)") & DERIV_MAX WRITE (N_PRT, "('AT ELEM =', I6, ', POINT = ', I2)") & TEST_IE, TEST_IP MAX_DERIVATIVE = DERIV_MAX ! FOR POSSIBLE PLOT SCALE IF ( U_PLT4 > 0 ) CLOSE (U_PLT4) !b IF ( U_PLT4 > 0 ) THEN !b CALL INQUIRE_ABOUT_UNIT (U_PLT4 ) !b PRINT *,'el_qp_xyz_fluxes CLOSED AT ELEMENT ', IE !b CLOSE (U_PLT4) !b END IF END IF ! LAST ELEMENT END IF ! ELEMENT OR BOUNDARY SEGMENT ! FLAG IF ANY SCP DATA WERE SAVED FOR AVERAGING CALL UPDATE_SCP_STATUS IF ( IE == N_ELEMS .AND. IO_1 == 0 ) PRINT *, & 'NOTE: SAVED ELEMENT FLUXES IN el_qp_xyz_fluxes.tmp' END SUBROUTINE LIST_ELEM_FLUXES SUBROUTINE LIST_ELEM_AND_EXACT_FLUXES (N_FILE, IE) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! LIST ELEMENT FLUXES AT QUADRATURE POINTS, ON N_FILE ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants Use Select_Source Use Elem_Type_Data ! for COORD (LT_N, N_SPACE), D (LT_FREE) IMPLICIT NONE ! ALWAYS USED INTEGER, INTENT(IN) :: N_FILE, IE ! Automatic Arrays !b REAL(DP) :: DGH (N_SPACE, LT_N), STRAIN (N_R_B + 2), & !b REAL(DP) :: STRAIN (N_R_B + 2), & !b XYZ (N_SPACE), E (N_R_B, N_R_B), STRESS (N_R_B + 2) REAL(DP) :: B (N_R_B, LT_FREE), FLUX (N_R_B), RMS !b INTEGER, SAVE :: TEST_IE, TEST_IP, J, N_IP, EOF, IO_1 REAL(DP), SAVE :: DERIV_MAX = 0.d0 ! VARIABLES: ! COORD = SPATIAL COORDINATES OF ELEMENT'S NODES ! D = NODAL PARAMETERS ASSOCIATED WITH AN ELEMENT ! E = ELEMENT CONSTITUTIVE MATRIX AT GAUSS POINT ! FLUX = EXACT FLUX COMPONENTS AT THE POINT ! G = GEOMETRIC INTERPOLATION FUNCTIONS ! H = SOLUTION INTERPOLATION FUNCTIONS ! LT_N = NUMBER OF NODES PER ELEMENT ! LT_FREE = NUMBER OF DEGREES OF FREEDOM PER ELEMENT ! N_ELEMS = TOTAL NUMBER OF ELEMENTS ! N_R_B = NUMBER OF ROWS IN B AND E MATRICES ! N_SPACE = DIMENSION OF SPACE ! N_FILE = UNIT FOR POST SOLUTION MATRICES STORAGE ! RMS = ROOT MEAN SQUARE VALUE ! STRAIN = STRAIN OR FLUX VECTOR ! XYZ = SPACE COORDINATES AT A POINT ! U_PLT4 = UNIT TO STORE PLOT DATA, IF > 0 !--> PRINT TITLES ON THE FIRST CALL IF ( IE == 1) THEN FLUX = 0.d0 REWIND (N_FILE) WRITE (N_PRT, 4) ; 4 FORMAT (/, & '*** FE AND EXACT FLUX COMPONENTS AT INTEGRATION POINTS ***',/, & '(SAVED FOR AVERAGING OR APPLICATION POST PROCESSING)' ) RECORD_NUMBER = 0 ! OPEN FLUX PLOT FILE IF ACTIVE IF ( U_PLT4 > 0 ) OPEN (U_PLT4, FILE='el_qp_xyz_fluxes.tmp', & ACTION='WRITE', STATUS='REPLACE', IOSTAT = IO_1) IF ( IO_1 > 0 ) THEN PRINT *, 'WARNING: OPEN FAILED FOR el_qp_xyz_fluxes.tmp' N_WARN = N_WARN + 1 END IF ENDIF ! IS THIS AN ELEMENT, OR A BOUNDARY SEGMENT ? IF ( IE <= N_ELEMS ) THEN ! ELEMENT RESULTS IF ( LIST_QP_FLUX ) THEN WRITE (N_PRT, 5) XYZ_SHORT (1:N_SPACE), FLUX_SHORT (1:N_R_B) !b WRITE (N_PRT, 6) XYZ_SHORT (1:N_SPACE), EXACT_SHORT (1:N_R_B) !b 5 FORMAT ('ELEMENT, PT, ', (9A10) ) 6 FORMAT (' PT, ', (9A10) ) END IF ! READ NUMBER OF INTEGRATION POINTS TO BE EMPLOYED READ (N_FILE, IOSTAT = EOF) N_IP IF ( EOF /= 0 ) THEN ! HIT END OF FILE PRINT *, 'LIST_ELEM_AND_EXACT_FLUXES EOF AT ELEMENT ', IE, EOF PRINT *, 'WAS EXPECTING: N_IP. WAS UNIT ', N_FILE, ' REWOUND ?' PRINT *, 'MAY NEED KEYWORD no_scp_ave IN CONTROL, OR MAY NEED' PRINT *, 'TO DELETE KEYWORD list_exact_flux, OR MAY NEED TO HAVE' PRINT *, 'STORE_FLUX_POINT_COUNT CALL IN my_el_sq_inc' CALL INQUIRE_ABOUT_UNIT (N_FILE) STOP 'ERROR, EOF IN LIST_ELEM_AND_EXACT_FLUXES' END IF !--> READ COORDS, CONSTITUTIVE, AND DERIVATIVE MATRIX DO J = 1, N_IP ! OVER ALL INTEGRATION POINTS READ (N_FILE, IOSTAT = EOF) XYZ, E, B IF ( EOF /= 0 ) THEN PRINT *, 'LIST_ELEM_AND_EXACT_FLUXES EOF AT ELEMENT ', IE, EOF PRINT *, 'WAS EXPECTING: XYZ, E, B AT POINT ', J, ' OF ', N_IP PRINT *, 'AT UNIT ', N_FILE PRINT *, 'MAY NEED KEYWORD no_scp_ave IN CONTROL, OR MAY NEED' PRINT *, 'TO DELETE KEYWORD list_exact_flux, OR MAY NEED' PRINT *, 'STORE_FLUX_POINT_DATA CALL IN my_el_sq_inc' CALL INQUIRE_ABOUT_UNIT (N_FILE) STOP 'ERROR, EOF IN LIST_ELEM_AND_EXACT_FLUXES' END IF ! CALCULATE DERIVATIVES, STRAIN = B * D STRAIN (1:N_R_B) = MATMUL (B, D) !b ! FLUX FROM CONSTITUTIVE DATA STRESS (1:N_R_B) = MATMUL (E, STRAIN (1:N_R_B)) if ( DEBUG_POST_EL .AND. EXAMPLE == 202 .AND. THIS_EL < 3 ) then print *, 'E ', THIS_EL; CALL RPRINT (E, N_R_B, N_R_B, 1) print *, 'STRAIN ' ; CALL RPRINT (STRAIN, 1, N_R_B, 1) print *, 'STRESS ' ; CALL RPRINT (STRESS, 1, N_R_B, 1) end if ! SCALAR MEASURE (VALID FOR VECTORS) RMS = SQRT ( SUM (STRESS (1:N_R_B)**2) ) ! GET EXACT FLUX VALUES CALL SELECT_EXACT_FLUX (XYZ, FLUX) IF ( DEBUG_B .AND. IE == 1 ) THEN PRINT *,'J ', J PRINT *,'E ' ; CALL RPRINT (E, N_R_B, N_R_B, 1) PRINT *,'B ' ; CALL RPRINT (B, N_R_B, LT_FREE, 1) END IF ! debug !--> PRINT COORDINATES AND FLUX AT THE POINT IF ( LIST_QP_FLUX ) THEN WRITE (N_PRT, '(I8, I3, 10(1PE11.3), (8X, 10(1PE11.3)))') & IE, J, XYZ, STRESS (1:N_R_B) WRITE (N_PRT, '(8X, I3, 10(1PE11.3), (8X, 10(1PE11.3)))') & J, XYZ, FLUX END IF IF ( RMS > DERIV_MAX ) THEN DERIV_MAX = RMS ! LARGEST VALUE TEST_IE = IE ; TEST_IP = J ! AT ELEM AND POINT END IF !--> STORE FLUX RESULTS TO BE PLOTTED LATER, IF USED IF ( U_PLT4 > 0 .AND. IO_1 == 0 ) & WRITE (U_PLT4, '( (10(1PE13.5)) )') XYZ, STRESS (1:N_R_B) ! SAVE COORDINATES & FLUX COMPONENTS FOR SCP FLUX AVERAGING IF ( U_SCPR > 0 ) THEN RECORD_NUMBER = RECORD_NUMBER + 1 SCP_RECORD_NUMBER (IE, J) = RECORD_NUMBER IF (RECORD_NUMBER == 1 .AND. DEBUG_SCP) & PRINT *,'Debug N_R_B, SCP_FIT ', N_R_B, SCP_FIT IF ( GRAD_BASE_ERROR ) THEN !b JEA 4/6/02 XXX WRITE (U_SCPR, REC = RECORD_NUMBER) XYZ, STRAIN (1:SCP_FIT) ELSE WRITE (U_SCPR, REC = RECORD_NUMBER) XYZ, STRESS (1:N_R_B) !b !b WRITE (U_SCPR, REC = RECORD_NUMBER) XYZ, STRESS (1:SCP_FIT) END IF ! GRAD VS FLUX END IF ! SCP RECOVERY END DO ! OVER INTEGRATION POINTS !b print *,'rec numb ', RECORD_NUMBER !b !b PRINT *, ' ' !b IF ( DEBUG_SCP ) CALL INQUIRE_ABOUT_SCP_UNIT !b !--> ARE FLUX CALCULATIONS COMPLETE FOR ALL ELEMENTS IF ( IE == N_ELEMS ) THEN ! LIST LARGEST FLUX WRITE (N_PRT, "('LARGEST FLUX RMS_VALUE = ', 1PE12.4)") & DERIV_MAX WRITE (N_PRT, "('AT ELEM =', I6, ', POINT = ', I2)") & TEST_IE, TEST_IP MAX_DERIVATIVE = DERIV_MAX ! FOR POSSIBLE PLOT SCALE IF ( U_PLT4 > 0 ) CLOSE (U_PLT4) END IF ! LAST ELEMENT END IF ! ELEMENT OR BOUNDARY SEGMENT ! FLAG IF ANY SCP DATA WERE SAVED FOR AVERAGING CALL UPDATE_SCP_STATUS IF ( IE == N_ELEMS .AND. IO_1 == 0 ) PRINT *, & 'NOTE: SAVED ELEMENT FLUXES (ONLY) IN el_qp_xyz_fluxes.tmp' END SUBROUTINE LIST_ELEM_AND_EXACT_FLUXES SUBROUTINE LIST_ELEM_CONVECTION_LOSS (N_FILE, IE, EACH) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! LIST CONVECTION HEAT LOSS IN EACH ELEMENT, ASSUMING ! CONSTANT CONVECTION COEFFICIENT AND EXTERNAL TEMPERATURE ! LOSS = integral ( h * H * (T - T_e)) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants Use Elem_Type_Data ! for D (LT_FREE) IMPLICIT NONE INTEGER, INTENT(IN) :: N_FILE, IE LOGICAL, INTENT(IN) :: EACH ! list each if true REAL(DP), SAVE :: VALUE, TOTAL, T_E REAL(DP) :: H_INTG (LT_FREE) ! integral h * H INTEGER :: EOF ! D = NODAL PARAMETERS ASSOCIATED WITH AN ELEMENT ! EACH = .TRUE. TO LIST EVERY VALUE, ELSE GIVE ONLY TOTAL ! H = SOLUTION INTERPOLATION FUNCTIONS ! H_INTG = INTEGRAL OF INTERPOLATION FUNCTIONS ! IE = ELEMENT, OR FACE NUMBER ! LT_FREE = NUMBER OF DEGREES OF FREEDOM PER ELEMENT ! N_FILE = UNIT FOR POST SOLUTION MATRICES STORAGE ! T_E = EXTERNAL TEMPERATURE FOR THE ELEMENT !--> PRINT TITLES ON THE FIRST CALL AND INITIALIZE IF ( IE == 1) THEN REWIND (N_FILE) TOTAL = 0.d0 IF ( EACH ) WRITE (N_PRT, 5) 5 FORMAT (/,'** CONVECTIVE HEAT LOSS RESULTS **',/, & 'ELEMENT LOSS') ENDIF !b print *, 'N_FILE ', N_FILE !b !b call LIST_SYSTEM_UNIT_NUMBERS !b IF ( IE <= N_ELEMS ) THEN ! ELEMENT RESULTS !--> GET ELEMENT EXTERNAL TEMPERATURE & CONVECTION INTEGRAL READ (N_FILE, IOSTAT = EOF) T_E, H_INTG IF ( EOF /= 0 ) THEN PRINT *, 'LIST_ELEM_CONVECTION_LOSS EOF AT ELEMENT ', IE PRINT *, 'WAS EXPECTING: T_E, H_INTG' STOP 'ERROR, EOF IN LIST_ELEM_CONVECTION_LOSS' END IF !--> CALCULATE ELEMENT CONTRIBUTION VALUE = DOT_PRODUCT (H_INTG, (D - T_E)) TOTAL = TOTAL + VALUE IF ( EACH ) WRITE (N_PRT, '(I7, 1PE18.6)') IE, VALUE IF (IE == N_ELEMS) WRITE (N_PRT, & "(/, 'TOTAL CONVECTIVE HEAT LOSS = ', 1PE16.6, /)") TOTAL END IF END SUBROUTINE LIST_ELEM_CONVECTION_LOSS SUBROUTINE LIST_REACTIONS_AT_ELEMS (NODES, DD) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! ELEMENT LEVEL REACTION RECOVERY ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants Use Elem_Type_Data Use Sys_Properties_Data ! for DELETED Use Interface_Header IMPLICIT NONE ! Material option for future uses. Add to arguments, etc. ! ALWAYS USED REAL(DP), INTENT(IN) :: DD (N_D_FRE) INTEGER, INTENT(IN) :: NODES (L_S_TOT, NOD_PER_EL) ! Automatic Arrays INTEGER :: IE, LT ! VARIABLES: ! C = ELEMENT INTERNAL SOURCE VECTOR ! D = NODAL PARAMETERS ASSOCIATED WITH AN ELEMENT ! DD = ARRAY OF SYSTEM DEGREES OF FREEDOM ! INDEX = SYSTEM DOF NOS ASSOCIATED WITH ELEMENT ! ELEM_NODES = THE NOD_PER_EL INCIDENCES OF THE ELEMENT ! NOD_PER_EL = NUMBER OF NODES PER ELEMENT ! N_D_FRE = TOTAL NUMBER OF SYSTEM DEGREES OF FREEDOM ! N_ELEMS = NUMBER OF ELEMENTS IN SYSTEM ! N_G_DOF = NUMBERS PARAMETERS PER NODE ! NODES = ELEMENT INCIDENCES OF ALL ELEMENTS ! L_REACT = UNIT FOR POST SOLUTION MATRICES STORAGE ! S = ELEMENT SQUARE MATRIX ! L_REACT MUST BE > 0 IF ( L_REACT == 0 ) THEN STOP 'NO L_REACT IN LIST_REACTIONS_AT_ELEMS' ELSE IF ( L_REACT > 0 ) REWIND (L_REACT) END IF LT = 1 ! INITIALIZE ELEMENT TYPES !--> LOOP OVER ELEMENTS DO IE = 1, N_ELEMS ! for elements !--> GET ELEMENT TYPE NUMBER IF ( N_L_TYPE > 1) LT = L_TYPE (IE) ! SAME AS LAST TYPE ? IF ( LT /= LAST_LT ) THEN ! this is a new type ! GET CONTROLS FOR THIS TYPE CALL SET_ELEM_TYPE_INFO (LT) END IF ! a new element type ! INITIALIZE AND CHECK IF ACTIVE IF ( DELETED (IE) ) CYCLE ! to skip deleted element !--> EXTRACT ELEMENT NODE NUMBERS ELEM_NODES = GET_ELEM_NODES (IE, LT_N, NODES) !9 !--> CALCULATE DEGREE OF FREEDOM NUMBERS INDEX = GET_ELEM_INDEX (LT_N, ELEM_NODES) !9 !--> EXTRACT NODAL PARAMETERS OF THE ELEMENT !b CALL EL_FRE (N_D_FRE, LT_FREE, D, DD, INDEX) !b CALL GET_ELEM_DOF (D, DD, INDEX) D = GET_ELEM_DOF (DD) !--> LIST THE REACTIONS, ELEMENT SOURCES, AND THEIR SUMS !b CALL GET_EL_REACTIONS (IE, LT_N, LT_FREE, N_G_DOF, & !b L_REACT, ELEM_NODES, D) IF ( LIST_LR > 0 .AND. L_REACT > 0 ) & CALL GET_EL_REACTIONS (IE, LT_N, LT_FREE, N_G_DOF, & !b L_REACT, ELEM_NODES, D) !b END DO ! over all elements and boundary segments !--> COMPLETED, FREE TYPE ARRAYS !b CALL DEALLOCATE_TYPE_INTERPOLATIONS !bb contains SUBROUTINE GET_EL_REACTIONS (IE, LT_N, LT_FREE, N_G_DOF, & N_FILE, ELEM_NODES, D) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! GET REACTIONS (FLUXES) AT AN ELEMENT'S NODES, ! R = S*D - C ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use Precision_Module IMPLICIT NONE INTEGER, INTENT(IN) :: IE, LT_N, LT_FREE, N_G_DOF INTEGER, INTENT(IN) :: N_FILE INTEGER, INTENT(IN) :: ELEM_NODES (LT_N) REAL(DP), INTENT(IN) :: D (LT_FREE) INTEGER :: IG, IN, IROW, EOF REAL(DP) :: S (LT_FREE, LT_FREE), C (LT_FREE) REAL(DP) :: ELEM_REACTS (LT_FREE) REAL(DP) :: SUM_REACT (N_G_DOF) REAL(DP) :: SUM_SOURCE (N_G_DOF) ! C = ELEMENT COLUMN MATRIX ! D = KNOWN SOLUTION CAUSING THE REACTIONS ! IE = CURRENT ELEMENT NUMBER ! ELEM_NODES = ELEMENT TOPOLOGY LIST ! ELEM_REACTS = ELEMENT REACTIONS VECTOR ! LT_N = NUMBER OF NODES PER ELEMENT ! LT_FREE = NUMBER OF ELEMENT DOF, N_G_DOF*N ! N_G_DOF = NUMBER OF GENERALIZED UNKNOWNS PER NODE ! N_FILE = UNIT TO HOLDING S & C FROM ELEMENT IE ! S = ELEMENT SQUARE MATRIX ! SUM_SOURCE = ELEMENT WORKSPACE RESULT ! SUM_REACT = NODAL WORKSPACE RESULT ! READ STIFFNESS, AND SOURCE FOR ELEM REACTIONS IF ( N_FILE < 1) STOP 'INVALID UNIT IN GET_EL_REACTIONS' IF ( IE == 1 ) WRITE (N_PRT, & "(/, '** ELEMENT REACTION, INTERNAL SOURCES AND SUMMATIONS **')") READ (N_FILE, IOSTAT = EOF) S, C IF ( EOF /= 0 ) THEN PRINT *, 'GET_EL_REACTIONS: EOF AT ELEMENT ', IE PRINT *, 'WAS EXPECTING: S, C' STOP 'ERROR, EOF IN GET_EL_REACTIONS' END IF WRITE (N_PRT, * ) ' ELEMENT ', IE WRITE (N_PRT, * ) ' NODE DOF REACTION ELEM_SOURCE SUMS' ! GIVE REACTIONS, R = S*D - C ELEM_REACTS = MATMUL (S, D) - C SUM_REACT = 0. ; SUM_SOURCE = 0. DO IN = 1, LT_N DO IG = 1, N_G_DOF IROW = N_G_DOF * (IN - 1) + IG SUM_SOURCE (IG) = SUM_SOURCE (IG) + C (IROW) SUM_REACT (IG) = SUM_REACT (IG) + ELEM_REACTS (IROW) WRITE (N_PRT, "( I5, I5, 1PE15.5, 1PE15.5 )") & ELEM_NODES (IN), IG, ELEM_REACTS (IROW), C (IROW) END DO END DO DO IG = 1, N_G_DOF PRINT "(' SUM:', I5, 3(1PE15.5))", IG, SUM_REACT (IG), & SUM_SOURCE (IG), (SUM_REACT (IG) + SUM_SOURCE (IG)) END DO END SUBROUTINE GET_EL_REACTIONS END SUBROUTINE LIST_REACTIONS_AT_ELEMS !b see library file_utility_lib.f for INQUIRE_ABOUT_UNIT !b see library file_utility_lib.f for INQUIRE_ABOUT_FILE SUBROUTINE LIST_SYSTEM_UNIT_NUMBERS Use System_Constants PRINT *, ' ' PRINT *, '** CURRENT SYSTEM UNIT NUMBER VALUES ** (0 IF UNUSED)' PRINT *, ' ' PRINT *, 'BINARY USER ON/OFF CONTROL:' PRINT *, 'USER APPLICATION OPTION ........ N_FILE1 = ', N_FILE1 PRINT *, 'USER APPLICATION OPTION ........ N_FILE2 = ', N_FILE2 PRINT *, 'USER APPLICATION OPTION ........ N_FILE3 = ', N_FILE3 PRINT *, 'USER APPLICATION OPTION ........ N_FILE4 = ', N_FILE4 PRINT *, 'USER APPLICATION OPTION ........ N_FILE5 = ', N_FILE5 PRINT *, ' ' PRINT *, 'FORMATTED I/O DEFAULTS:' PRINT *, 'STANDARD INPUT ................. N_CRD = ', N_CRD PRINT *, 'STANDARD OUTPUT ................ N_PRT = ', N_PRT PRINT *, 'DEBUG OUTPUT ................... N_BUG = ', N_BUG PRINT *, ' ' PRINT *, 'FORMATTED RESULTS (AND MATLAB PLOTS):' PRINT *, 'RESULT dof_history_###.tmp ..... U_DOF = ', U_DOF PRINT *, 'RESULT node_history_###.tmp .... U_NODE = ', U_NODE PRINT *, 'RESULT msh_bc_xyz.tmp .......... U_PLT1 = ', U_PLT1 PRINT *, 'RESULT msh_typ_nodes.tmp ....... U_PLT2 = ', U_PLT2 PRINT *, 'RESULT node_results.tmp ........ U_PLT3 = ', U_PLT3 PRINT *, 'RESULT el_qp_xyz_fluxes.tmp .... U_PLT4 = ', U_PLT4 PRINT *, 'RESULT scp_node_ave_fluxes.tmp . U_PLT5 = ', U_PLT5 PRINT *, 'RESULT el_error_est.tmp ........ U_PLT6 = ', U_PLT6 PRINT *, 'RESULT pt_ave_error_est.tmp .... U_PLT7 = ', U_PLT7 PRINT *, 'RESULT pt_ave_grad_flux.tmp .... U2_PLT1 = ', U2_PLT1 PRINT *, 'EXACT pt_ex_grad_flux.tmp .... U2_PLT2 = ', U2_PLT2 PRINT *, 'EXACT exact_node_solution.tmp .. U_PLT8 = ', U_PLT8 PRINT *, 'EXACT exact_node_flux.tmp ...... U_PLT9 = ', U_PLT9 PRINT *, 'Note: See CALL TURN_OFF_PLOTS, and keywords sav_*' PRINT *, ' ' PRINT *, 'BINARY INTERNAL DEFAULTS:' PRINT *, 'SYSTEM REACTION STORAGE ........ N_REACT = ', N_REACT PRINT *, 'ELEMENTS REACTION STORAGE ...... L_REACT = ', L_REACT PRINT *, 'SUPER_CONVERGENT_PATCH ......... U_SCPR = ', U_SCPR PRINT *, 'SCP 2ND DERIVATIVES ............ U_2_DER = ', U_2_DER PRINT *, 'SCP 3RD DERIVATIVES ............ U_3_DER = ', U_3_DER PRINT *, 'Note: See CALL TURN_ON_SCP, CALL TURN_OFF_SCP, no_scp_ave' PRINT *, 'ELEMENT FLUX DATA .............. U_FLUX = ', U_FLUX PRINT *, 'ELEMENT INTEGRATIONS ........... U_INTGR = ', U_INTGR PRINT *, 'Note: CALL INQUIRE_ABOUT_UNIT (unit) for current status' PRINT *, 'Note: CALL LIST_OPEN_IO_UNITS for current open status' PRINT *, ' ' PRINT *, 'FORMATED JOB LOG, AND HISTORIES' PRINT *, 'CURRENT JOB LOG ................ U_LOG = ', U_LOG PRINT *, 'HISTORY OF RECENT JOB LOGS ..... U_HIST = ', U_HIST PRINT *, 'keyword_input.echo ............ECHO_UNIT = ', ECHO_UNIT END SUBROUTINE LIST_SYSTEM_UNIT_NUMBERS SUBROUTINE READ_STATUS (IO_STATUS, ITEM, SUB_ITEM) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! IF ERROR, LIST IO_STATUS, ITEM AND SUB_ITEM NUMBERS ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants ! For N_WARN, LINE_WARN IMPLICIT NONE INTEGER, INTENT (IN) :: IO_STATUS, ITEM, SUB_ITEM !b INTEGER, PARAMETER :: LIMIT = 5 INTEGER, SAVE :: COUNT = 0 SELECT CASE ( IO_STATUS ) CASE (:-1) ! END OF FILE CALL CITE_INPUT_LINE ! List current input position PRINT *, 'INPUT ERROR: END OF FILE AT ITEM ', ITEM IF (SUB_ITEM > 0 ) PRINT *, 'AND SUB_ITEM ', SUB_ITEM STOP 'INPUT ERROR: END OF FILE' CASE (1:) ! PROBABLE USER FORMAT ERROR CALL CITE_INPUT_LINE ! List current input position PRINT *, 'WARNING: USER FORMAT ERROR LIKELY AT ITEM ', ITEM IF (SUB_ITEM > 0 ) PRINT *, 'AND SUB_ITEM ', SUB_ITEM N_WARN = N_WARN + 1 COUNT = COUNT + 1 CASE DEFAULT ; RETURN ! NO ERROR END SELECT !b IF ( COUNT == LIMIT ) STOP 'READ_STATUS: MAX WARNINGS REACHED' IF ( COUNT == LINE_WARN ) THEN PRINT *,'ERROR: MAX WARNINGS FOR 1 LINE REACHED, ENTER keyword' PRINT *,'LINE_WARN WITH A VALUE > ', LINE_WARN STOP 'READ_STATUS: MAX WARNINGS REACHED' END IF ! avoid infinite loops END SUBROUTINE READ_STATUS SUBROUTINE LIST_2ND_DERIV_ORDER ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! List the storage order used for second derivatives in MODEL ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants PRINT *,'NOTE, EXACT AND FE 2ND DERIVATIVES ARE STORED HERE AS:' SELECT CASE (N_SPACE) CASE (1) ; PRINT *,'EXACT: U,xx FE: U,xx' CASE (2) ; PRINT *,'EXACT 3: U,xx U,yx=U,xy U,yy' PRINT *,'FE 4: U,xx U,yx U,xy U,yy' CASE (3) PRINT *,'EXACT 6: U,xx U,yx=U,xy U,zx=U,xz U,yy U,zy=U,yz U,zz' PRINT *,'FE 9: U,xx U,yx U,zx U,xy U,yy U,zy U,xz U,yz U,zz' CASE DEFAULT PRINT *,'ARE NOT DEFINED IN THIS DIMENSION SPACE' END SELECT END SUBROUTINE LIST_2ND_DERIV_ORDER SUBROUTINE GET_STEP_RESULTS (DD) ! * * * * * * * * * * * * * * * * * * * * * * * * * * ! GET NODAL RESULTS DATA FOR RESTART ! Requires Global RESTART_STEP >= 0 ! * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants Use Interface_Header IMPLICIT NONE CHARACTER (LEN=3) :: T_STRING ! for file names CHARACTER (LEN=20) :: F_NAME ! file name REAL(DP), INTENT(OUT) :: DD (N_D_FRE) INTEGER :: I, K, IO_1, U_TMP, GET_NEXT_IO_UNIT ! Automatic Arrays INTEGER :: INDEX (N_G_DOF) ! DD = CALCULATED NODAL PARAMETERS ! MAX_NP = NUMBER OF NODES IN SYSTEM ! N_D_FRE = NUMBER OF DOF IN THE SYSTEM ! N_G_DOF = NUMBER OF PARAMETERS (DOF) PER NODE ! U_TMP = UNIT TO OUTPUT NODE RESULTS FOR PLOTS U_TMP = GET_NEXT_IO_UNIT () T_STRING = TO3STR (RESTART_STEP) F_NAME = 'node_results_' // T_STRING // '.tmp' OPEN (U_TMP, FILE=F_NAME, ACTION='READ', & STATUS='OLD', IOSTAT = IO_1) IF ( IO_1 > 0 ) THEN PRINT *, 'WARNING: D==0, RESTART OPEN FAILED. ', F_NAME N_WARN = N_WARN + 1 ; DD = 0.d0 ; RETURN END IF ! OPEN FAILURE DO I = 1, MAX_NP INDEX = GET_INDEX_AT_PT (I) READ (U_TMP, *) (DD (INDEX (K)), K = 1, N_G_DOF) !b 5 FORMAT ( (10(1X, 1PE12.5)) ) END DO CLOSE (U_TMP) ; PRINT *, ' ' PRINT *, 'NOTE: GOT RESTART SOLUTION FROM ', F_NAME END SUBROUTINE GET_STEP_RESULTS SUBROUTINE SAVE_STEP_RESULTS (DD) ! * * * * * * * * * * * * * * * * * * * * * * * * * * ! SAVE NODAL RESULTS DATA FOR PLOTS ! * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants Use Interface_Header IMPLICIT NONE CHARACTER (LEN=3) :: T_STRING ! for file names CHARACTER (LEN=20) :: F_NAME ! file name REAL(DP), INTENT(IN) :: DD (N_D_FRE) INTEGER :: I, K, IO_1, U_TMP, GET_NEXT_IO_UNIT ! Automatic Arrays INTEGER :: INDEX (N_G_DOF) ! DD = CALCULATED NODAL PARAMETERS ! MAX_NP = NUMBER OF NODES IN SYSTEM ! N_D_FRE = NUMBER OF DOF IN THE SYSTEM ! N_G_DOF = NUMBER OF PARAMETERS (DOF) PER NODE ! U_TMP = UNIT TO OUTPUT NODE RESULTS FOR PLOTS U_TMP = GET_NEXT_IO_UNIT () T_STRING = TO3STR (THIS_STEP) F_NAME = 'node_results_' // T_STRING // '.tmp' OPEN (U_TMP, FILE=F_NAME, ACTION='WRITE', & STATUS='REPLACE', IOSTAT = IO_1) IF ( IO_1 > 0 ) THEN PRINT *, 'WARNING: OPEN FAILED. ', F_NAME N_WARN = N_WARN + 1 ; RETURN END IF ! OPEN FAILURE DO I = 1, MAX_NP INDEX = GET_INDEX_AT_PT (I) WRITE (U_TMP, 5) (DD (INDEX (K)), K = 1, N_G_DOF) 5 FORMAT ( (10(1X, 1PE12.5)) ) END DO CLOSE (U_TMP) PRINT *, 'NOTE: SOLUTION RESULTS SAVED TO ', F_NAME PRINT *, ' ' END SUBROUTINE SAVE_STEP_RESULTS SUBROUTINE SAVE_STEP_RESULTS_WITH_EXACT (DD, X) ! * * * * * * * * * * * * * * * * * * * * * * * * * * ! SAVE NODAL AND EXACT RESULTS DATA FOR PLOTS ! * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants Use Select_Source Use Interface_Header IMPLICIT NONE CHARACTER (LEN=3) :: T_STRING ! for file names CHARACTER (LEN=18) :: F_NAME ! file name REAL(DP), INTENT(IN) :: DD (N_D_FRE) REAL(DP), INTENT(IN) :: X (MAX_NP, N_SPACE) INTEGER :: I, K, IO_1, U_TMP, GET_NEXT_IO_UNIT ! Automatic Arrays INTEGER :: INDEX (N_G_DOF) REAL(DP) :: EXACT_SOL (N_G_DOF) ! DD = CALCULATED NODAL PARAMETERS ! EXACT_SOL = USER SUPPLIED SOLUTION AT A POINT ! MAX_NP = NUMBER OF NODES IN SYSTEM ! N_D_FRE = NUMBER OF DOF IN THE SYSTEM ! N_G_DOF = NUMBER OF PARAMETERS (DOF) PER NODE ! U_PLT3 = UNIT TO OUTPUT NODE COORDINATES FOR PLOTS ! X = COORDINATES OF ALL SYSTEM NODES !b IF ( U_PLT3 < 1 ) RETURN ! output turned off U_TMP = GET_NEXT_IO_UNIT () T_STRING = TO3STR (THIS_STEP) F_NAME = 'node_exact_' // T_STRING // '.tmp' OPEN (U_TMP , FILE=F_NAME, ACTION='WRITE', & STATUS='REPLACE', IOSTAT = IO_1) IF ( IO_1 > 0 ) THEN PRINT *, 'WARNING: OPEN FAILED. ', F_NAME N_WARN = N_WARN + 1 ; RETURN END IF ! OPEN FAILURE DO I = 1, MAX_NP !b CALL EXACT_SOLUTION (X(I, :), EXACT_SOL) CALL SELECT_EXACT_SOLUTION (X(I, :), EXACT_SOL) !b INDEX = GET_INDEX_AT_PT (I) WRITE (U_TMP , 5) (DD (INDEX (K)), K = 1, N_G_DOF), & EXACT_SOL 5 FORMAT ( (10(1X, 1PE12.5)) ) END DO CLOSE (U_TMP ) PRINT *, 'NOTE: FE & EXACT RESULTS SAVED TO ', F_NAME END SUBROUTINE SAVE_STEP_RESULTS_WITH_EXACT SUBROUTINE SAVE_STEP_EXACT_FLUXES_AT_NODES (X) ! * * * * * * * * * * * * * * * * * * * * * * * * * * ! SAVE EXACT FLUXES TO A FILE ! * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants Use Select_Source IMPLICIT NONE CHARACTER (LEN=3) :: T_STRING ! for file names CHARACTER (LEN=23) :: F_NAME ! file name REAL(DP), INTENT(IN) :: X (MAX_NP, N_SPACE) REAL(DP) :: FLUX (N_R_B) INTEGER :: IP, IO_1, U_TMP, GET_NEXT_IO_UNIT ! FLUX = EXACT FLUX COMPONENTS AT A POINT ! MAX_NP = NUMBER OF SYSTEM NODES ! X = COORDINATES OF ALL SYSTEM NODES ! U_PLT9 = UNIT TO OUTPUT NODAL GRADIENTS, IF U_SCPR>0 IF ( .NOT. USE_EXACT_FLUX ) THEN PRINT *, 'WARNING, SAVE_STEP_EXACT_FLUXES_AT_NODES: ', & 'NO EXACT SOLUTION USED' N_WARN = N_WARN + 1 ; RETURN END IF !b IF ( U_PLT9 < 1 ) RETURN ! output turned off U_TMP = GET_NEXT_IO_UNIT () T_STRING = TO3STR (THIS_STEP) F_NAME = 'exact_node_flux_' // T_STRING // '.tmp' OPEN (U_TMP, FILE=F_NAME, ACTION='WRITE', & STATUS='REPLACE', IOSTAT = IO_1) IF ( IO_1 > 0 ) THEN PRINT *, 'OPEN FAILED FOR ', F_NAME N_WARN = N_WARN + 1 ; RETURN END IF DO IP = 1, MAX_NP !b CALL EXACT_FLUX (X (IP, :), FLUX) CALL SELECT_EXACT_FLUX (X (IP, :), FLUX) !b WRITE (U_TMP, 5) FLUX 5 FORMAT ( (8(1PE13.5)) ) END DO ! FOR ALL NODES CLOSE (U_TMP) PRINT *, 'NOTE: EXACT NODAL FLUXES SAVED TO ', F_NAME END SUBROUTINE SAVE_STEP_EXACT_FLUXES_AT_NODES SUBROUTINE SAVE_STEP_EXACT_NODE_SOL (X) ! * * * * * * * * * * * * * * * * * * * * * * * * * * ! SAVE NODAL AND EXACT RESULTS DATA FOR PLOTS ! * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants Use Select_Source Use Interface_Header IMPLICIT NONE CHARACTER (LEN=3) :: T_STRING ! for file names CHARACTER (LEN=22) :: F_NAME ! file name REAL(DP), INTENT(IN) :: X (MAX_NP, N_SPACE) INTEGER :: I, IO_1 ! Automatic Arrays REAL(DP) :: EXACT_SOL (N_G_DOF) ! EXACT_SOL = USER SUPPLIED SOLUTION AT A POINT ! MAX_NP = NUMBER OF NODES IN SYSTEM ! N_D_FRE = NUMBER OF DOF IN THE SYSTEM ! N_G_DOF = NUMBER OF PARAMETERS (DOF) PER NODE ! U_PLT8 = UNIT TO OUTPUT NODE SOLUTION FOR PLOTS ! X = COORDINATES OF ALL SYSTEM NODES IF ( U_PLT8 < 1 ) RETURN ! output turned off T_STRING = TO3STR (THIS_STEP) F_NAME = 'exact_node_sol_' // T_STRING // '.tmp' OPEN (U_PLT8, FILE=F_NAME, ACTION='WRITE', & STATUS='REPLACE', IOSTAT = IO_1) IF ( IO_1 > 0 ) THEN PRINT *,'WARNING: OPEN FAILED ', F_NAME N_WARN = N_WARN + 1 ; RETURN END IF ! OPEN FAILURE DO I = 1, MAX_NP !b CALL EXACT_SOLUTION (X (I, :), EXACT_SOL) CALL SELECT_EXACT_SOLUTION (X (I, :), EXACT_SOL) !b WRITE (U_PLT8, 5) EXACT_SOL 5 FORMAT ( (10(1X, 1PE12.5)) ) END DO CLOSE (U_PLT8) PRINT *, 'NOTE: EXACT NODAL VALUES SAVED TO ', F_NAME END SUBROUTINE SAVE_STEP_EXACT_NODE_SOL