! File: model_lib.f ! copyright 2005, J. E. Akin, all rights reserved. ! --------------------------------------------------------------- ! Alphabetical listings follow ! --------------------------------------------------------------- SUBROUTINE B_MATRIX_EQUAL_DGH (DGH, B) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * ! B == DGH FOR MOST SCALAR PROBLEMS, SO JUST COPY IT ! * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants Use Elem_Type_Data ! for LT_N, LT_FREE IMPLICIT NONE REAL (DP), INTENT(IN) :: DGH (N_SPACE, LT_N) REAL (DP), INTENT(OUT) :: B (N_R_B, LT_FREE) IF ( N_SPACE == N_R_B .AND. LT_N == LT_FREE ) THEN B (1:N_SPACE, 1:LT_N) = DGH ELSE ! INVALID CALL PRINT *,'ERROR: B_MATRIX_EQUAL_DGH, INCOMPATABLE MATRIX SIZES' PRINT *,'EXPECTED SIZE: ', N_R_B, ' BY ', LT_FREE PRINT *,'RECEIVED SIZE: ', N_SPACE, ' BY ', LT_N STOP 'INCOMPATABLE MATRIX SIZES, B_MATRIX_EQUAL_DGH' END IF END SUBROUTINE B_MATRIX_EQUAL_DGH !SUBROUTINE BOUNDARY_FLUX (FLUX, COORD_SEG, C_SEG, S_SEG, I_OPT, & ! NSQP, NSPARM, ISEG) !! * * * * * * * * * * * * * * * * * * * * * * * * * * !! PROBLEM DEPENDENT BOUNDARY FLUX CONTRIBUTIONS !! >>> Outdated version, inactive in F90 version !! * * * * * * * * * * * * * * * * * * * * * * * * * * !Use System_Constants !Use Sys_Properties_Data !Use Elem_Type_Data ! for: ! ! LT_FREE, LT_GEOM, LT_N, LT_PARM, LT_QP, ELEM_NODES (LT_N), & ! ! COORD (LT_N, N_SPACE), GEOMETRY (LT_GEOM, N_SPACE), & ! ! C (LT_FREE), D (LT_FREE), S (LT_FREE, LT_FREE), & ! ! EL_M (LT_FREE, LT_FREE), DIAG (LT_FREE), & ! ! DLG (LT_PARM, LT_GEOM), DLG_QP (LT_PARM, LT_GEOM, LT_QP) & ! ! DLH (LT_PARM, LT_N), DLH_QP (LT_PARM, LT_N, LT_QP), & ! ! DLV (LT_PARM, LT_FREE), DLV_QP (LT_PARM, LT_FREE, LT_QP), & ! ! G (LT_GEOM), G_QP (LT_GEOM, LT_QP), & ! ! H (LT_N), H_QP (LT_N, LT_QP), & ! ! V (LT_FREE), V_QP (LT_FREE, LT_QP), & ! ! PT (LT_PARM, LT_QP), WT (LT_QP) !Use Interface_Header ! for functions ! IMPLICIT NONE !! ALWAYS USED ! INTEGER, INTENT(IN) :: NSQP, NSPARM, ISEG ! REAL(DP), INTENT(IN) :: COORD_SEG (L_B_N, N_SPACE) ! REAL(DP), INTENT(IN) :: FLUX (L_B_N, N_G_FLUX) ! REAL(DP), INTENT(OUT) :: C_SEG (N_D_FLUX), S_SEG (N_D_FLUX, N_D_FLUX) ! INTEGER, INTENT(OUT) :: I_OPT !! Automatic Arrays ! !b REAL(DP) :: XYZ (N_SPACE), AJ (LT_PARM, LT_PARM), & ! !b DGH (N_SPACE, LT_N), DGV (N_SPACE, LT_FREE), & ! !b AJ_INV (LT_PARM, LT_PARM), GPT (LT_QP), GWT (LT_QP) ! !! AJ = JACOBIAN !! AJ_INV = JACOBIAN INVERSE !! C_SEG = BOUNDARY FLUX COLUMN MATRIX CONTRIBUTIONS !! COORD_SEG = SPATIAL COORDINATES OF SEGMENT NODES !! DGH = GLOBAL DERIVATIVES SCALAR INTERPOLATION FUNCTIONS !! DGV = GLOBAL DERIVATIVES VECTOR INTERPOLATION FUNCTIONS !! DLG = LOCAL DERIVATIVES GEOMETRIC INTERPOLATION !! DLH = LOCAL DERIVATIVES SCALAR INTERPOLATION FUNCTIONS !! DLV = LOCAL DERIVATIVES VECTOR INTERPOLATION FUNCTIONS !! FLT_SP = REAL PROPERTIES ON THE SEGMENTS !! FLUX = SPECIFIED BOUNDARY FLUX COMPONENTS !! G = GEOMETRIC INTERPOLATION FUNCTIONS !! H = SOLUTION INTERPOLATION FUNCTIONS !! I_OPT = PROBLEM MATRIX REQUIREMENTS, MUST BE SET. !! = 1, CALCULATE C ONLY !! = 2, CALCULATE S ONLY !! = 3, CALCULATE BOTH C AND S !! ISEG = SEGMENT NUMBER !! L_B_N = NO. OF NODES ON AN ELEMENT BOUNDARY SEGMENT !! L_HOMO = 1 IF SEGMENT PROPERTIES ARE HOMOGENEOUS !! NOD_PER_EL = NUMBER OF SOLUTION NODES, NOD_PER_EL = L_B_N USUALLY !! N_BS_FIX = NUMBER OF INTEGER PROPERTIES PER SEGMENT !! N_BS_FLO = NUMBER OF REAL PROPERTIES PER SEGMENT !! SP_FIX = INTEGER PROPERTIES ON THE SEGMENTS !! N_D_FLUX = L_B_N*N_G_DOF = MAXIMUM NUMBER OF FLUX CONTRIBUTIONS !! N_G_DOF = NUMBER OF PARAMETERS PER NODE POINT !! N_PARM = PARAMETRIC GEOMETRY NODES, = N_SPACE USUALLY !! N_QP = NUMBER OF QUADRATURE POINTS !! N_F_SEG = MAX NUMBER OF SEGMENTS !! N_SPACE = DIMENSION OF SOLUTION SPACE !! PT = QUADRATURE COORDINATES !! S_SEG = BOUNDARY FLUX SQUARE MATRIX !! V = VECTOR INTERPOLATION FUNCTIONS !! WT = QUADRATURE WEIGHTS !! XYZ = SPACE COORDINATES AT A POINT !! .................................................... !! ** BOUNDARY_FLUX PROBLEM DEPENDENT STATEMENTS FOLLOW ** !! .................................................... !! I_OPT = PROBLEM MATRIX REQUIREMENTS, MUST BE SET. ! S_SEG = 0.d0 ; C_SEG = 0.d0 ! I_OPT = 0 ! or REPLACE with: include 'my_flux_inc' !! use print to supress compiler warnings !if ( I_BUG > 0 ) print *, NSQP, NSPARM, ISEG !b !if ( I_BUG > 0 ) print *, COORD_SEG !b !if ( I_BUG > 0 ) print *, FLUX !b !! include 'my_flux_inc' !END SUBROUTINE BOUNDARY_FLUX SUBROUTINE CALC_PROP_AT_PT (H, PRT_L_PT, VALUES) ! * * * * * * * * * * * * * * * * * * * * * * * * * * ! CALCULATE N_NP_FLO PROPERTIES AT A LOCAL PT USING ! ELEMENT'S NODAL PROPERTIES, PRT_L_PT, AND THE N ! INTERPOLATION FUNCTIONS, H, AT THE POINT ! * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants IMPLICIT NONE REAL(DP), INTENT(IN) :: H (NOD_PER_EL) REAL(DP), INTENT(IN) :: PRT_L_PT (NOD_PER_EL, N_NP_FLO) REAL(DP), INTENT(OUT) :: VALUES (N_NP_FLO) ! NOD_PER_EL = NUMBER OF NODES PER ELEMENT ! N_NP_FLO = NO. OF FLOATING POINT NODAL PROPERTIES ! H = C^0 INTERPOLATION FUNCTIONS FOR AN ELEMENT ! PRT_L_PT = FLOATING PT PROPS OF ELEMENT'S NODES ! VALUES = LOCAL VALUES OF PROPERTIES IF ( N_NP_FLO < 1 ) STOP 'pt_real = 0, CALC_PROP_AT_PT' VALUES = MATMUL (H, PRT_L_PT) END SUBROUTINE CALC_PROP_AT_PT SUBROUTINE CHANGE_IN_SOLUTION (DD, DD_OLD, OLD_NORM, NEW_NORM, & DIFF_NORM, RATIO) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! CALCULATE THE MEAN CHANGE IN NODAL PARAMETERS FROM ! THE LAST ITERATION ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants IMPLICIT NONE REAL(DP), INTENT(IN) :: DD (N_D_FRE), DD_OLD(N_D_FRE) REAL(DP), INTENT(OUT) :: OLD_NORM, NEW_NORM REAL(DP), INTENT(OUT) :: DIFF_NORM, RATIO REAL(DP), PARAMETER :: EPS = EPSILON (1.d0) ! * CHANGE SHOULD BE CALLED BEFORE CORRECT_ITERATION * ! RATIO = DIFF_NORM/NEW_NORM ! DIFF_NORM = SQRT(SUM OF (DD(I) - DD_OLD(I))**2) ! OLD_NORM = SQRT(SUM OF DD_OLD(I)**2) ! NEW_NORM = SQRT(SUM OF DD(I)**2) ! DD_OLD = NODAL PARAMETER LIST FROM LAST ITERATION ! DD = NODAL PARAMETERS FROM CURRENT ITERATION ! N_D_FRE = TOTAL NO OF DEGREES OF FREEDOM IN SYS NEW_NORM = SQRT ( SUM (DD**2) ) OLD_NORM = SQRT ( SUM (DD_OLD**2) ) DIFF_NORM = SQRT ( SUM ((DD - DD_OLD) **2) ) RATIO = DIFF_NORM / (NEW_NORM + EPS) WRITE (N_PRT, 5) OLD_NORM, NEW_NORM, DIFF_NORM, RATIO 5 FORMAT ( & '*** NODAL DOF FOR CURRENT AND PREVIOUS ITERATIONS ***',/,& 'ROOT MEAN SQ OF PREVIOUS VALUES . . . ',1PE13.5,/, & 'ROOT MEAN SQ OF CURRENT VALUES . . . ',1PE13.5,/, & 'ROOT MEAN SQ OF DIFFERENCES . . . . . ',1PE13.5,/, & 'RATIO OF LAST TWO QUANTITIES . . . . ',1PE13.5,/) END SUBROUTINE CHANGE_IN_SOLUTION SUBROUTINE CHECK_EL_SHAPE ! * * * * * * * * * * * * * * * * * * * * * * * * * * * ! CHECK L_SHAPE DATA FOR FREQUENT USER ERRORS ! * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants IMPLICIT NONE LOGICAL :: WARN ! L_B_N = NUMBER OF BOUNDARY NODES ON ELEM ! L_SHAPE = ELEMENT SHAPE CODE ! NOD_PER_EL = NUMBER OF NODES PER ELEMENT ! N_SPACE = DIMENSION OF ELEMENT SPACE ! >>> Note the CONTINUITY constant is now available here too WARN = .FALSE. ! FORCE THE INPUT VALUE IF NEGATIVE IF ( L_SHAPE < 0 ) THEN L_SHAPE = IABS (L_SHAPE) PRINT *, 'WARNING, L_SHAPE SIGN CORRECTED' N_WARN = N_WARN + 1 ! INCREMENT WARNING RETURN END IF ! LINE IF ( NOD_PER_EL == 2 .AND. L_SHAPE /= 1 ) THEN L_SHAPE = 1 ; WARN = .TRUE. END IF ! TRIANGLE IF ( N_SPACE > 1 .AND. NOD_PER_EL == 3 .AND. L_SHAPE /= 2 ) THEN L_SHAPE = 2 ; WARN = .TRUE. END IF ! QUADRILATERAL IF ( N_SPACE > 1 .AND. NOD_PER_EL == 4 .AND. L_SHAPE /= 3 ) THEN L_SHAPE = 3 ; WARN = .TRUE. END IF ! TRI OR WEDGE IF ( NOD_PER_EL == 6 .AND. L_SHAPE /= 2 ) THEN IF ( N_SPACE == 2 ) THEN L_SHAPE = 2 ; WARN = .TRUE. END IF END IF ! QUAD OR HEX IF ( NOD_PER_EL == 8 .AND. L_SHAPE /= 3 ) THEN IF ( N_SPACE /= 3 ) THEN L_SHAPE = 3 ; WARN = .TRUE. END IF END IF IF ( WARN ) THEN WRITE (N_PRT, 1) L_SHAPE 1 FORMAT ('********************************',/, & '* WARNING L_SHAPE CHANGED TO', I2, ' *', /, & '********************************') N_WARN = N_WARN + 1 ! INCREMENT WARNING END IF WARN = .FALSE. IF ( L_B_N > 0 ) THEN ! CHECK EDGE_FACE_VOLUME DATA (INCOMPLETE, NEEDS L_SHAPE ?) !bb xx IF (NOD_PER_EL == 2 .AND. L_B_N /= 1) WARN = .TRUE. IF (NOD_PER_EL == 3 .AND. L_B_N /= 2) WARN = .TRUE. IF (NOD_PER_EL == 4 .AND. L_B_N /= 2) WARN = .TRUE. IF (NOD_PER_EL == 6 .AND. L_B_N /= 3) WARN = .TRUE. IF (NOD_PER_EL == 8 .AND. L_B_N > 4) WARN = .TRUE. IF ( WARN ) THEN PRINT *, 'WARNING: CHECK_EL_SHAPE, el_segment VALUE ', L_B_N PRINT *, 'UNEXPECTED FOR el_nodes VALUE ', NOD_PER_EL N_WARN = N_WARN + 1 ! INCREMENT WARNING END IF END IF WARN = .FALSE. IF ( M_B_N > 0 ) THEN ! CHECK EDGE_FACE_VOLUME DATA (INCOMPLETE, NEEDS L_SHAPE ?) !bb xx IF (NOD_PER_EL == 2 .AND. M_B_N /= 1) WARN = .TRUE. IF (NOD_PER_EL == 3 .AND. M_B_N /= 2) WARN = .TRUE. IF (NOD_PER_EL == 4 .AND. M_B_N /= 2) WARN = .TRUE. IF (NOD_PER_EL == 6 .AND. M_B_N /= 3) WARN = .TRUE. IF (NOD_PER_EL == 8 .AND. M_B_N > 4) WARN = .TRUE. IF ( WARN ) THEN PRINT *, 'WARNING: CHECK_EL_SHAPE, MIXED NODES VALUE ', M_B_N PRINT *, 'UNEXPECTED FOR el_nodes VALUE ', NOD_PER_EL N_WARN = N_WARN + 1 ! INCREMENT WARNING END IF END IF END SUBROUTINE CHECK_EL_SHAPE SUBROUTINE CONSTRAINT_COUNT (N_RES_EQ, I_BC) ! * * * * * * * * * * * * * * * * * * * * * * * * * * ! CALCULATE NUMBER OF CONSTRAINT FLAGS OF EACH TYPE ! * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants IMPLICIT NONE INTEGER, INTENT(IN) :: I_BC (MAX_NP) INTEGER, INTENT(INOUT) :: N_RES_EQ (MAX_TYP) INTEGER :: I, ITEST, J, K ! Automatic Arrays !b INTEGER :: KODES (N_G_DOF) INTEGER :: KODES (N_G_DOF), USED (N_G_DOF) ! MAX_NP = TOTAL NUMBER OF SYSTEM NODES ! N_G_DOF = NO. OF PARAMETERS (DOF) PER NODE ! I_BC = NODAL POINT BOUNDARY RESTRAINT INDICATOR ! KODES = LIST OF RESTRAINT INDICATORS AT A NODE ! N_RES_EQ = NUMBER OF CONSTR EQS ON EXIT, N_RES_EQ ! MAX_TYP = MAX NO OF DIFFERENT CONSTRAINT TYPES ! MAX_ACT = ACTIVE NO OF TYPES ! USED = DATA CHECKING WORK ARRAY USED = 0 ! INITIALIZATION N_RES_EQ = 0 ! INITIALIZATION DO I = 1, MAX_NP ! DOES NODE I HAVE A NODAL PARAMETER CONSTRAINT ITEST = IABS (I_BC (I) ) IF ( ITEST > 0 ) THEN ! EXTRACT PARAMETER CODES CALL PT_CODE (I, ITEST, KODES) DO J = 1, N_G_DOF K = KODES (J) ! UPDATE CONSTRAINT COUNTERS !b IF ( K > 0) N_RES_EQ (K) = N_RES_EQ (K) + 1 IF ( K > 0 ) THEN N_RES_EQ (K) = N_RES_EQ (K) + 1 USED (J) = 1 END IF ! USED END DO END IF END DO ! CONVERT TO EQUATION COUNTERS N_CEQ = 0 MAX_ACT = 0 WRITE (N_PRT, 5) 5 FORMAT ( /,'*** NODAL PARAMETER CONSTRAINT LIST ***', /,& 'TYPE EQUATIONS') DO I = 1, MAX_TYP K = N_RES_EQ (I) IF ( K > 0) MAX_ACT = I IF ( ( (K / I) *I) < K ) THEN WRITE (N_BUG, * ) 'INVALID DATA FOR CONSTRAINT TYPE', I N_WARN = N_WARN + 1 END IF N_RES_EQ (I) = N_RES_EQ (I) / I IF ( N_RES_EQ (I) > 0 ) WRITE (N_PRT, 5020) I, N_RES_EQ (I) 5020 FORMAT ( I4, I10 ) N_CEQ = N_CEQ + N_RES_EQ (I) END DO ! VERIFY ALL DOF USED !b DO J = 1, N_G_DOF !b IF ( USED (J) <= 0 ) THEN !b PRINT *,'WARNING: PARAMETER ', J, ' HAS NO BOUNDARY CONDITION' !b N_WARN = N_WARN + 1 !b END IF ! UNUSED IN BC !b END DO ! VALIDATION !b !b IF ( ABS(N_CEQ) > N_D_FRE ) STOP 'impossible constraints' IF ( ABS(N_CEQ) > N_D_FRE ) THEN PRINT *,'ERROR: CONSTRAINT_COUNT EXCEEDS NUMBER OF EQUATIONS' STOP 'ERROR: CONSTRAINT_COUNT EXCEEDS NUMBER OF EQUATIONS' END IF ! IMPOSSIBLE DATA END SUBROUTINE CONSTRAINT_COUNT FUNCTION COPY_DGH_INTO_B_MATRIX (DGH) RESULT (B) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * ! B == DGH FOR MOST SCALAR PROBLEMS, SO JUST COPY IT ! * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants Use Elem_Type_Data ! for LT_N, LT_FREE IMPLICIT NONE REAL (DP), INTENT(IN) :: DGH (N_SPACE, LT_N) REAL (DP) :: B (N_R_B, LT_FREE) IF ( N_SPACE == N_R_B .AND. LT_N == LT_FREE ) THEN B (1:N_SPACE, 1:LT_N) = DGH ELSE ! INVALID CALL STOP 'INCOMPATABLE MATRIX SIZES, COPY_DGH_INTO_B_MATRIX' END IF END FUNCTION COPY_DGH_INTO_B_MATRIX SUBROUTINE PUT_DGH_INTO_B_MATRIX (DGH, B) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * ! B == DGH FOR MOST SCALAR PROBLEMS, SO JUST COPY IT ! * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants Use Elem_Type_Data ! for LT_N, LT_FREE IMPLICIT NONE REAL (DP), INTENT(IN) :: DGH (N_SPACE, LT_N) REAL (DP) :: B (N_R_B, LT_FREE) IF ( N_SPACE == N_R_B .AND. LT_N == LT_FREE ) THEN B (1:N_SPACE, 1:LT_N) = DGH ELSE ! INVALID CALL STOP 'INCOMPATABLE MATRIX SIZES, PUT_DGH_INTO_B_MATRIX' END IF END SUBROUTINE PUT_DGH_INTO_B_MATRIX SUBROUTINE CORRECT_ITERATION (DD, DD_OLD) ! * * * * * * * * * * * * * * * * * * * * * * * * * * ! CALCULATE NEW STARTING VALUES FOR NEXT ITERATION ! * * * * * * * * * * * * * * * * * * * * * * * * * * ! OVER RELAXATION METHOD Use System_Constants IMPLICIT NONE REAL(DP), INTENT(IN) :: DD (N_D_FRE) REAL(DP), INTENT(INOUT) :: DD_OLD (N_D_FRE) REAL(DP) :: OMEGA (N_G_DOF) INTEGER :: I, J, K ! DD = CALCULATED DOF FROM LAST ITERATION ! DD_OLD = DOF TO BE USED TO START NEXT ITERATION ! N_D_FRE = TOTAL NO OF SYS DEGREES OF FREEDOM ! N_G_DOF = NUMBER OF PARAMETERS (DOF) PER NODE ! OMEGA = RELAXATION FACTOR FOR EACH NODAL DOF !b OMEGA = 1.25_DP ! TYPICAL VALUE IF ( SPACE_TIME .AND. (.NOT. S_T_CONTINUOUS) ) THEN IF ( MAX_NP /= (MAX_NP/2)*2 ) THEN ! not space time PRINT *,'WARNING, CORRECT_ITERATION: Odd max node number' N_WARN = N_WARN + 1 END IF ! move t_(n+1) answers to t_(n) for weak BC arrays I = N_D_FRE / 2 DD_OLD (1:I) = DD_OLD (I+1:I+N_D_FRE/2) ELSE ! standard iteration OMEGA = RELAXATION ! TYPICAL VALUE I = 0 DO J = 1, MAX_NP DO K = 1, N_G_DOF I = I + 1 DD_OLD (I) = DD_OLD (I) + OMEGA (K) * (DD (I) - DD_OLD (I) ) END DO END DO END IF END SUBROUTINE CORRECT_ITERATION FUNCTION COUNT_EL_TYPE_NODES (LT) RESULT (COUNT) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! GET DATA BASED ON ELEMENT TYPE (SEE CONTROL_INPUT ) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants Use Elem_Type_Data ! for LT_DATA (MAX_DAT, N_L_TYPE) IMPLICIT NONE INTEGER, INTENT(IN) :: LT INTEGER :: COUNT ! LT = ELEMENT TYPE NUMBER ! N_L_TYPE = NUMBER OF ELEMENT TYPES ! LT_N = NUMBER OF NODES PER ELEMENT LT_N = LT_DATA (1, LT) COUNT = LT_N END FUNCTION COUNT_EL_TYPE_NODES !SUBROUTINE DEGPAR (I_PT, J_PARM, N_G_DOF, INDEX) !! * * * * * * * * * * * * * * * * * * * * * * * !! DETERMINE THE DEGREE OF FREEDOM NUMBER !! OF NODAL PARAMETER J_PARM AT NODE POINT I_PT !! * * * * * * * * * * * * * * * * * * * * * * * !Use Precision_Module ! IMPLICIT NONE ! INTEGER, INTENT(IN) :: I_PT, J_PARM, N_G_DOF ! INTEGER, INTENT(OUT) :: INDEX ! !! I_PT = A SPECIFIC NODE NUMBER !! J_PARM = THE PARAMETER AT THAT NODE !! N_G_DOF = NUMBER OF PARAMETERS PER NODE ! ! !b IF ( DOF_VARY ) STOP 'DEGPAR: INVALID FOR DOF_VARY' ! ! INDEX = N_G_DOF * (I_PT - 1) + J_PARM !END SUBROUTINE DEGPAR SUBROUTINE ELEM_COORD (NOD_ON_EL, L_SPACE, X, COORD, ELEM_NODES) ! * * * * * * * * * * * * * * * * * * * * * * * * * * ! DETERMINE COORDINATES OF NODES ON ELEMENT ! * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants IMPLICIT NONE INTEGER, INTENT(IN) :: NOD_ON_EL, L_SPACE REAL(DP), INTENT(IN) :: X (MAX_NP, L_SPACE) INTEGER, INTENT(IN) :: ELEM_NODES (NOD_ON_EL) REAL(DP), INTENT(OUT) :: COORD (NOD_ON_EL, L_SPACE) INTEGER :: I ! MAX_NP = NUMBER OF NODES IN SYSTEM ! L_SPACE = DIMENSION OF SPACE ! NOD_ON_EL = NUMBER OF NODES PER ELEMENT ! X = COORDINATES OF SYSTEM NODES ! COORD = COORDINATES OF ELEMENT NODES ! ELEM_NODES = NOD_ON_EL INCIDENCES OF ELEMENT DO I = 1, NOD_ON_EL ! ALLOW FOR OMITTED NODES IF ( ELEM_NODES (I) > 0) & COORD (I, 1:L_SPACE) = X (ELEM_NODES (I), 1:L_SPACE) END DO END SUBROUTINE ELEM_COORD FUNCTION GET_ELEM_COORD (X, ELEM_NODES) RESULT (COORD) ! * * * * * * * * * * * * * * * * * * * * * * * * * * ! DETERMINE COORDINATES OF NODES ON ELEMENT ! * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants IMPLICIT NONE REAL(DP), INTENT(IN) :: X (MAX_NP, N_SPACE) INTEGER, INTENT(IN) :: ELEM_NODES (NOD_PER_EL) REAL(DP) :: COORD (NOD_PER_EL, N_SPACE) INTEGER :: I ! MAX_NP = NUMBER OF NODES IN SYSTEM ! N_SPACE = DIMENSION OF SPACE ! NOD_PER_EL = NUMBER OF NODES PER ELEMENT ! X = COORDINATES OF SYSTEM NODES ! COORD = COORDINATES OF ELEMENT NODES ! ELEM_NODES = NOD_PER_EL INCIDENCES OF ELEMENT DO I = 1, NOD_PER_EL ! ALLOW FOR OMITTED NODES IF ( ELEM_NODES (I) > 0) & COORD (I, 1:N_SPACE) = X (ELEM_NODES (I), 1:N_SPACE) END DO END FUNCTION GET_ELEM_COORD SUBROUTINE ELEMENT_NODES (L_ID, LT_N, NODES, ELEM_NODES) ! * * * * * * * * * * * * * * * * * * * * * * * * ! EXTRACT NODES ASSOCIATED WITH ELEMENT L_ID ! * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants INTEGER, INTENT(IN) :: NODES (L_S_TOT, NOD_PER_EL), L_ID, LT_N INTEGER, INTENT(OUT) :: ELEM_NODES (LT_N) ! N_ELEMS = NUMBER OF ELEMENTS IN SYSTEM ! NOD_PER_EL = MAXIMUN NUMBER OF NODES PER ELEMENT ! LT_N = NUMBER OF NODES PER ELEMENT ! L_ID = ELEMENT NUMBER ! NODES = NODAL INCIDENCES OF ALL ELEMENTS ! ELEM_NODES = THE NOD_PER_EL NODAL INCIDENCES OF THE ELEMENT IF ( LT_N > NOD_PER_EL ) STOP 'DATA ERROR ELEMENT_NODES' ELEM_NODES (1:LT_N) = NODES (L_ID, 1:LT_N) END SUBROUTINE ELEMENT_NODES SUBROUTINE EL_FRE (N_D_FRE, N_EL_FRE, D, DD, INDEX) ! * * * * * * * * * * * * * * * * * * * * * * * * * * ! EXTRACT ELEMENT DEGREES OF FREEDOM FROM SYSTEM DOF ! * * * * * * * * * * * * * * * * * * * * * * * * * * Use Precision_Module IMPLICIT NONE INTEGER, INTENT(IN) :: N_D_FRE, N_EL_FRE, INDEX (N_EL_FRE) REAL(DP), INTENT(IN) :: DD (N_D_FRE) REAL(DP), INTENT(OUT) :: D (N_EL_FRE) ! D = NODAL PARAMETERS ASSOCIATED WITH ELEMENT ! DD = SYSTEM ARRAY OF NODAL PARAMETERS ! INDEX = ARRAY OF SYSTEM DEGREE OF FREEDOM NUMBERS ! N_EL_FRE = NUMBER OF DEGREES OF FREEDOM PER ELEMENT ! N_D_FRE = TOTAL NUMBER OF SYSTEM DEGREES OF FREEDOM WHERE ( INDEX > 0 ) D = DD (INDEX) ELSEWHERE D = 0.d0 END WHERE END SUBROUTINE EL_FRE FUNCTION GET_ELEM_DOF (DD) RESULT (D_LT) ! * * * * * * * * * * * * * * * * * * * * * * * * * * ! EXTRACT ELEMENT DEGREES OF FREEDOM FROM SYSTEM DOF ! * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants ! for N_D_FRE Use Elem_Type_Data ! for LT_FREE, INDEX (LT_FREE) IMPLICIT NONE REAL(DP), INTENT(IN) :: DD (N_D_FRE) REAL(DP) :: D_LT (LT_FREE) ! D_LT = NODAL PARAMETERS ASSOCIATED WITH ELEMENT ! DD = SYSTEM ARRAY OF NODAL PARAMETERS ! INDEX = ARRAY OF SYSTEM DEGREE OF FREEDOM NUMBERS ! LT_FREE = NUMBER OF DEGREES OF FREEDOM FOR ELEMENT TYPE ! N_D_FRE = TOTAL NUMBER OF SYSTEM DEGREES OF FREEDOM WHERE ( INDEX > 0 ) D_LT = DD (INDEX) ELSEWHERE D_LT = 0.d0 END WHERE END FUNCTION GET_ELEM_DOF SUBROUTINE GET_TYPE_DOF (LT_D, DD, LT_INDEX) ! * * * * * * * * * * * * * * * * * * * * * * * * * * ! GET ELEMENT TYPE DEGREES OF FREEDOM FROM SYSTEM DOF ! * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants ! for N_D_FREE Use Elem_Type_Data ! for LT_FREE IMPLICIT NONE INTEGER, INTENT(IN) :: LT_INDEX (LT_FREE) REAL(DP), INTENT(IN) :: DD (N_D_FRE) REAL(DP), INTENT(OUT) :: LT_D (LT_FREE) ! LT_D = NODAL PARAMETERS ASSOCIATED WITH ELEMENT ! DD = SYSTEM ARRAY OF NODAL PARAMETERS ! LT_INDEX = ARRAY OF SYSTEM DEGREE OF FREEDOM NUMBERS ! LT_FREE = NUMBER OF DEGREES OF FREEDOM FOR ELEMENT TYPE ! N_D_FRE = TOTAL NUMBER OF SYSTEM DEGREES OF FREEDOM WHERE ( LT_INDEX > 0 ) LT_D = DD (LT_INDEX) ELSEWHERE LT_D = 0.d0 END WHERE END SUBROUTINE GET_TYPE_DOF SUBROUTINE ELEM_PROPERTIES (L_ID, L_PROP, EL_PROP) ! * * * * * * * * * * * * * * * * * * * * * * * * * * ! EXTRACT PROPERTIES OF ELEMENT, L_ID, FROM TOTAL ! PROPERTIES ARRAYS >>> now inactive, but correct ! * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants Use Sys_Properties_Data IMPLICIT NONE INTEGER, INTENT(IN) :: L_ID REAL(DP), INTENT(OUT) :: EL_PROP (N_LP_FLO) INTEGER, INTENT(OUT) :: L_PROP (N_LP_FIX) INTEGER :: I ! L_ID = ELEMENT NUMBER ! LP_FIX = SYSTEM ARRAY OF FIXED PT ELEM PROPERTIES ! L_PROP = ELEM FIXED PT PROPERTIES ARRAY ! FLT_EL = SYS ARRAY OF FLOATING PT NODAL PROP ! EL_PROP = ELEM FLOATING PT PROPERTIES ARRAY ! L_HOMO = 1, IF PROPERTIES ARE SAME IN ALL ELEMENTS ! N_LP_FIX = NUMBER OF INTEGER ELEMENT PROPERTIES ! N_LP_FLO = NUMBER OF REAL ELEMENT PROPERTIES I = L_ID IF ( L_HOMO == 1 ) I = 1 ! FLOATING POINT PROPERTIES EL_PROP (1:N_LP_FLO) = FLT_EL (I, 1:N_LP_FLO) ! FIXED POINT PROPERTIES L_PROP (1:N_LP_FIX) = LP_FIX (I, 1:N_LP_FIX) END SUBROUTINE ELEM_PROPERTIES SUBROUTINE GENERATE_EL_MATRICES_NEW (IE, X, DD_OLD) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! GENERATE ELEMENT MATRICES AND POST SOLUTION DATA ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants Use Elem_Type_Data ! for PT (LT_PARM, LT_QP), WT (LT_QP), ! G (LT_GEOM, LT_QP), H (LT_N ), ! DLG (LT_PARM, LT_GEOM, LT_QP), ! DLH (LT_PARM, LT_N , LT_QP), ! EL_M (LT_FREE, LT_FREE), DIAG_M (LT_FREE), ! S (LT_FREE, LT_FREE), C (LT_FREE), ! ELEM_NODES (LT_N), ! COORD (LT_N, N_SPACE), ! GEOMETRY (LT_GEOM, N_SPACE) ! INDEX (LT_FREE), D (LT_FREE) Use Sys_Properties_Data ! for FLT_MISC, MISC_FIX, DELETED Use Select_Source Use Interface_Header IMPLICIT NONE INTEGER, INTENT(IN) :: IE ! SYSTEM DATA REAL(DP), INTENT(IN) :: X (MAX_NP, N_SPACE), DD_OLD (N_D_FRE) ! Automatic Arrays REAL(DP) :: E (N_R_B, N_R_B) REAL(DP) :: H_INTG (LT_N ) REAL(DP) :: PRT_L_PT (LT_N, N_NP_FLO) REAL(DP) :: PRT_MAT (MISC_FL) REAL(DP) :: FLUX (L_B_N, N_G_FLUX) !b INTEGER :: L_PT_PROP (LT_N, N_NP_FIX) ! VARIABLES: ! COORD = SPATIAL COORDINATES OF ELEMENT'S NODES ! D = NODAL PARAMETERS ASSOCIATED WITH AN ELEMENT ! DD_OLD = SYSTEM NODAL PARAMETERS FROM LAST ITERATION ! DGH = GLOBAL DERIVATIVES INTERPOLATION FUNCTIONS ! DLG = LOCAL DERIVATIVES GEOMETRIC INTERPOLATION ! DLH = LOCAL DERIVATIVES INTERPOLATION FUNCTIONS ! E = CONSTITUTIVE MATRIX ! FLT_EL = REAL PROPERTIES OF ELEMENTS ! FLT_MISC = MISCELLANEOUS REAL PROPERTIES OF SYSTEM ! FLT_NP = REAL PROPERTIES OF SYSTEM NODES ! FLT_SP = REAL PROPERTIES OF SYSTEM SEGMENTS ! G = GEOMETRIC INTERPOLATION FUNCTIONS ! GEOMETRY = COORDINATES OF ELEM GEOMETRIC NODES ! H = SOLUTION INTERPOLATION FUNCTIONS ! H_INTG = INTEGRAL OF INTERPOLATION FUNCTIONS ! IE = CURRENT ELEMENT NUMBER ! INDEX = SYSTEM DOF NUMBERS ASSOCIATED WITH ELEMENT ! L_HOMO = 1, IF ELEMENT PROPERTIES ARE HOMOGENEOUS ! ELEM_NODES = THE NOD_PER_EL INCIDENCES OF THE ELEMENT ! LP_FIX = SYSTEM ARRAY OF INTEGER ELEM PROPERTIES ! L_PT_PROP = INTEGER PROPERTIES AT EACH ELEMENT NODE ! LP_TEST > 0, IF ELEMENT PROPERTIES HAVE BEEN DEFINED ! LT = ELEMENT TYPE NUMBER (IF USED) ! LT_USER = OPTIONAL USER CODE ASSOCIATED WITH TYPE ! MAX_NP = NUMBER OF SYSTEM NODES ! MISC_FIX = MISCELLANEOUS INTEGER SYSTEM PROPERTIES ! LT_N = NUMBER OF NODES PER ELEMENT ! N_D_FRE = TOTAL NUMBER OF SYSTEM DEGREES OF FREEDOM ! N_ELEMS = NUMBER OF ELEMENTS ! LT_FREE = NUMBER OF DEGREES OF FREEDOM PER ELEMENT ! LT_GEOM = NUMBER OF GEOMETRY NODES ! N_HOMO = 1, IF NODAL PROPERTIES ARE HOMOGENEOUS ! N_ITER = NO. OF ITERATIONS TO BE RUN (USUALLY 1) ! N_MAT = MATERIAL TYPE NUMBER ! LT_PARM = DIMENSION OF PARAMETRIC SPACE ! NP_FIX = INTEGER PROPERTIES AT ALL NODES ! LT_QP = NUMBER OF QUADRATURE POINTS ! N_R_B = NUMBER OF ROWS IN B AND E MATRICES ! N_SPACE = DIMENSION OF SPACE ! N_FILE1 = UNIT FOR POST SOLUTION MATRICES STORAGE ! N_FILE2 = OPTIONAL UNIT FOR USER (USED WHEN > 0) ! N_FILE# = OPTIONAL UNITS FOR USER (USED WHEN > 0) ! NUL_COL = 1, IF ELEMENT COLUMN MATRIX IS ALWAYS ZERO ! PRT_L_PT = REAL PROPERTIES AT ELEMENT NODES ! PRT_MAT = REAL ELEM PROPERTIES BASED ON MATERIAL NUMBER ! PT = QUADRATURE COORDINATES ! S = ELEMENT SQUARE MATRIX ! SP_FIX = INTEGER PROPERTIES AT ALL SEGMENTS ! WT = QUADRATURE WEIGHTS ! X = COORDINATES OF SYSTEM NODES IF ( I_BUG > 0 .AND. (N_ITER > 0 .OR. TRANSIENT .OR. DYNAMIC)) THEN PRINT *,'Enter GENERATE_EL_MATRICES_NEW with DD_OLD' !b CALL RPRINT (DD_OLD, 1, N_D_FRE, 1) END IF ! GET ELEMENT GROUP: STANDARD, MIXED BC, SEGMENT, USER DEFINED CALL SET_TO_MIXED_OR_SEG_OR_STD (IE) ! For IS_STD, IS_MIX, IS_SEG !--> EXTRACT ANALYSIS AND GEOMETRIC NODAL COORDINATES CALL ELEM_COORD (LT_N, N_SPACE, X, COORD, ELEM_NODES) CALL ELEM_COORD (LT_GEOM, N_SPACE, X, GEOMETRY, ELEM_NODES(1:LT_GEOM)) !--> EXTRACT NODAL PARAMETERS FROM LAST ITERATION (IF ANY) IF ( N_ITER > 1 .OR. TRANSIENT ) D = GET_ELEM_DOF (DD_OLD) !b call at(693); print *,'DD_OLD ', DD_OLD ; print *, 'D ', THIS_EL, D !--> EXTRACT NODAL POINT PROPERTIES (IF ANY) IF ( N_NP_FLO > 0) CALL ELEM_NODE_PROPERTIES (LT_N, PRT_L_PT, & L_PT_PROP, ELEM_NODES) !b call at(640) !--> EXTRACT REAL MATERIAL PROPERTIES (IF ANY) IF ( N_MAT > 0) CALL MATERIAL_PROP (N_MAT, PRT_MAT) !--> GENERATE ELEMENT SQUARE AND COLUMN MATRICES E = GET_REAL_IDENTITY (N_R_B) ! DEFAULT TO IDENTITY MATRIX H_INTG = 0.d0 ! ZERO INTEGRAL IF (IS_ELEMENT) THEN ! Get sq and maybe column matrices CALL SELECT_ELEM_SQ_MATRIX (E, H_INTG, PRT_L_PT, PRT_MAT, & L_PT_PROP, IE, X) !b for space frame IF ( DEBUG_EL_COL .AND. IE <= 3 ) CALL RPRINT &!b (C, 1, LT_FREE, 1) !b IF ( NUL_COL == 0 ) THEN ! Need column matrix CALL SELECT_ELEM_COL_MATRIX (E, H_INTG, PRT_L_PT, & PRT_MAT, L_PT_PROP, IE) END IF ! Non-zero column matrix IF ( DEBUG_EL_COL .AND. IE <= 3 ) CALL RPRINT &!b (C, 1, LT_FREE, 1) !b END IF ! Standard element IF (IS_MIXED_BC) CALL SELECT_MIXED_SQ_MATRIX (E, H_INTG, PRT_L_PT,& PRT_MAT, L_PT_PROP, IE) IF (IS_FLUX_SEG) THEN ! Get and use boundary flux input !--> EXTRACT FLUX IF BOUNDARY SEGMENT ELEMENT IF ( IS_SEG > 0 ) FLUX (:, :) = S_FLUX(:, :, IS_SEG) !b make call CALL SELECT_SEG_COL_MATRIX (E, H_INTG, PRT_L_PT, PRT_MAT, & L_PT_PROP, IE, FLUX) END IF ! Boundary flux source !--> STORE DATA FOR POST SOLUTION CALCULATIONS (IF ANY) IF ( N_FILE1 > 0 .OR. N_FILE2 > 0 ) & CALL SELECT_ELEM_POST_DATA (E, H_INTG, PRT_L_PT, PRT_MAT, & L_PT_PROP, IE, FLUX) ! NOTE: SYSTEM PROPERTIES UPDATE COULD BE DONE HERE END SUBROUTINE GENERATE_EL_MATRICES_NEW SUBROUTINE GENERATE_EL_MATRICES (IE, X, DD_OLD) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! GENERATE ELEMENT MATRICES AND POST SOLUTION DATA ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants Use Elem_Type_Data ! for PT (LT_PARM, LT_QP), WT (LT_QP), ! G (LT_GEOM, LT_QP), H (LT_N ), ! DLG (LT_PARM, LT_GEOM, LT_QP), ! DLH (LT_PARM, LT_N , LT_QP), ! S (LT_FREE, LT_FREE), C (LT_FREE), ! ELEM_NODES (LT_N), ! COORD (LT_N, N_SPACE), ! GEOMETRY (LT_GEOM, N_SPACE) ! INDEX (LT_FREE), D (LT_FREE) Use Sys_Properties_Data ! for FLT_MISC, MISC_FIX, DELETED Use Select_Source Use Interface_Header IMPLICIT NONE INTEGER, INTENT(IN) :: IE ! SYSTEM DATA REAL(DP), INTENT(IN) :: X (MAX_NP, N_SPACE), DD_OLD (N_D_FRE) ! Automatic Arrays REAL(DP) :: E (N_R_B, N_R_B) REAL(DP) :: H_INTG (LT_N ) !b REAL(DP) :: PRT_L_PT (LT_N, N_NP_FLO), EL_PROP (N_LP_FLO) REAL(DP) :: PRT_L_PT (LT_N, N_NP_FLO) REAL(DP) :: PRT_MAT (MISC_FL) REAL(DP) :: FLUX (L_B_N, N_G_FLUX) !b !b INTEGER :: L_PROP (N_LP_FIX), L_PT_PROP (LT_N, N_NP_FIX) INTEGER :: L_PT_PROP (LT_N, N_NP_FIX) !b INTEGER :: I_SEG !b ! VARIABLES: ! COORD = SPATIAL COORDINATES OF ELEMENT'S NODES ! D = NODAL PARAMETERS ASSOCIATED WITH AN ELEMENT ! DD_OLD = SYSTEM NODAL PARAMETERS FROM LAST ITERATION ! DGH = GLOBAL DERIVATIVES INTERPOLATION FUNCTIONS ! DLG = LOCAL DERIVATIVES GEOMETRIC INTERPOLATION ! DLH = LOCAL DERIVATIVES INTERPOLATION FUNCTIONS ! E = CONSTITUTIVE MATRIX ! FLT_EL = REAL PROPERTIES OF ELEMENTS ! FLT_MISC = MISCELLANEOUS REAL PROPERTIES OF SYSTEM ! FLT_NP = REAL PROPERTIES OF SYSTEM NODES ! FLT_SP = REAL PROPERTIES OF SYSTEM SEGMENTS ! G = GEOMETRIC INTERPOLATION FUNCTIONS ! GEOMETRY = COORDINATES OF ELEM GEOMETRIC NODES ! H = SOLUTION INTERPOLATION FUNCTIONS ! H_INTG = INTEGRAL OF INTERPOLATION FUNCTIONS ! IE = CURRENT ELEMENT NUMBER ! INDEX = SYSTEM DOF NUMBERS ASSOCIATED WITH ELEMENT ! L_HOMO = 1, IF ELEMENT PROPERTIES ARE HOMOGENEOUS ! ELEM_NODES = THE NOD_PER_EL INCIDENCES OF THE ELEMENT ! LP_FIX = SYSTEM ARRAY OF INTEGER ELEM PROPERTIES ! L_PT_PROP = INTEGER PROPERTIES AT EACH ELEMENT NODE ! LP_TEST > 0, IF ELEMENT PROPERTIES HAVE BEEN DEFINED ! LT = ELEMENT TYPE NUMBER (IF USED) ! LT_USER = OPTIONAL USER CODE ASSOCIATED WITH TYPE ! MAX_NP = NUMBER OF SYSTEM NODES ! MISC_FIX = MISCELLANEOUS INTEGER SYSTEM PROPERTIES ! LT_N = NUMBER OF NODES PER ELEMENT ! N_D_FRE = TOTAL NUMBER OF SYSTEM DEGREES OF FREEDOM ! N_ELEMS = NUMBER OF ELEMENTS ! LT_FREE = NUMBER OF DEGREES OF FREEDOM PER ELEMENT ! LT_GEOM = NUMBER OF GEOMETRY NODES ! N_HOMO = 1, IF NODAL PROPERTIES ARE HOMOGENEOUS ! N_ITER = NO. OF ITERATIONS TO BE RUN (USUALLY 1) ! N_MAT = MATERIAL TYPE NUMBER ! LT_PARM = DIMENSION OF PARAMETRIC SPACE ! NP_FIX = INTEGER PROPERTIES AT ALL NODES ! LT_QP = NUMBER OF QUADRATURE POINTS ! N_R_B = NUMBER OF ROWS IN B AND E MATRICES ! N_SPACE = DIMENSION OF SPACE ! N_FILE1 = UNIT FOR POST SOLUTION MATRICES STORAGE ! N_FILE2 = OPTIONAL UNIT FOR USER (USED WHEN > 0) ! N_FILE# = OPTIONAL UNITS FOR USER (USED WHEN > 0) ! NUL_COL = 1, IF ELEMENT COLUMN MATRIX IS ALWAYS ZERO ! PRT_L_PT = REAL PROPERTIES AT ELEMENT NODES ! PRT_MAT = REAL ELEM PROPERTIES BASED ON MATERIAL NUMBER ! PT = QUADRATURE COORDINATES ! S = ELEMENT SQUARE MATRIX ! SP_FIX = INTEGER PROPERTIES AT ALL SEGMENTS ! WT = QUADRATURE WEIGHTS ! X = COORDINATES OF SYSTEM NODES IF ( DEBUG_ASSEMBLY ) THEN PRINT *,'Enter GENERATE_EL_MATRICES with DD_OLD' CALL RPRINT (DD_OLD, 1, N_D_FRE, 1) END IF ! GET ELEMENT GROUP: STANDARD, MIXED BC, SEGMENT, USER DEFINED CALL SET_TO_MIXED_OR_SEG_OR_STD (IE) ! For IS_STD, IS_MIX, IS_SEG !--> EXTRACT ANALYSIS AND GEOMETRIC NODAL COORDINATES CALL ELEM_COORD (LT_N, N_SPACE, X, COORD, ELEM_NODES) CALL ELEM_COORD (LT_GEOM, N_SPACE, X, GEOMETRY, ELEM_NODES(1:LT_GEOM)) !--> EXTRACT NODAL PARAMETERS FROM LAST ITERATION (IF ANY) IF ( TRANSIENT .OR. N_ITER > 1) D = GET_ELEM_DOF (DD_OLD) !b call at(829); print *,'DD_OLD ', DD_OLD ; print *, 'D ', THIS_EL, D !--> EXTRACT NODAL POINT PROPERTIES (IF ANY) IF ( N_NP_FLO > 0) CALL ELEM_NODE_PROPERTIES (LT_N, PRT_L_PT, & L_PT_PROP, ELEM_NODES) !--> EXTRACT REAL MATERIAL PROPERTIES (IF ANY) IF ( N_MAT > 0) CALL MATERIAL_PROP (N_MAT, PRT_MAT) !--> GENERATE ELEMENT SQUARE AND COLUMN MATRICES E = GET_REAL_IDENTITY (N_R_B) ! DEFAULT TO IDENTITY MATRIX H_INTG = 0.d0 ! ZERO INTEGRAL IF (IS_ELEMENT ) THEN ! Get sq and maybe column matrices CALL ELEM_SQ_MATRIX (E, H_INTG, PRT_L_PT, PRT_MAT, & L_PT_PROP, IE, X) !b for space frame !b L_PT_PROP, IE) IF ( NUL_COL == 0 ) CALL ELEM_COL_MATRIX (E, H_INTG, PRT_L_PT, & PRT_MAT, L_PT_PROP, IE) END IF ! Standard element !b IF (IS_MIXED_BC) CALL MIXED_SQ_MATRIX (E, H_INTG, PRT_L_PT, & IF (IS_MIXED_BC) CALL SELECT_MIXED_SQ_MATRIX (E, H_INTG, PRT_L_PT,& PRT_MAT, L_PT_PROP, IE) IF (IS_FLUX_SEG) THEN !--> EXTRACT FLUX IF BOUNDARY SEGMENT ELEMENT IF ( IS_SEG > 0 ) FLUX (:, :) = S_FLUX(:, :, IS_SEG) CALL SEG_COL_MATRIX (E, H_INTG, PRT_L_PT, PRT_MAT, L_PT_PROP, & IE, FLUX) END IF !--> STORE DATA FOR POST SOLUTION CALCULATIONS (IF ANY) IF ( N_FILE1 > 0 .OR. N_FILE2 > 0 ) & CALL ELEM_POST_DATA (E, H_INTG, PRT_L_PT, PRT_MAT, & L_PT_PROP, IE, FLUX) ! NOTE: SYSTEM PROPERTIES UPDATE COULD BE DONE HERE END SUBROUTINE GENERATE_EL_MATRICES FUNCTION GET_PT_DOF_INDEX (I_PT, J_PARM) RESULT (INDEX) ! * * * * * * * * * * * * * * * * * * * * * * * ! DETERMINE THE DEGREE OF FREEDOM NUMBER ! OF NODAL PARAMETER J_PARM AT NODE POINT I_PT ! * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants ! for N_G_DOF IMPLICIT NONE INTEGER, INTENT(IN) :: I_PT, J_PARM INTEGER :: INDEX ! DOF_PT = NUMBER OF DOF AT EACH NODE: (0:MAX_NP) ! DOF_PT_SUM = TOTAL NUMBER OF DOF, INCLUDING THIS NODE: (-1:MAX_NP) ! I_PT = A SPECIFIC NODE NUMBER ! J_PARM = THE PARAMETER AT THAT NODE ! N_G_DOF = NUMBER OF PARAMETERS PER NODE IF ( DOF_VARY ) THEN ! VARIABLE DOF PER NODE INDEX = DOF_PT_SUM (I_PT-1) + J_PARM IF ( J_PARM > DOF_PT (I_PT) ) THEN PRINT *,'WARNING: INVALID DOF ', J_PARM, ' AT NODE ', I_PT N_WARN = N_WARN + 1 INDEX = 0 END IF ELSE ! ALL NODES HAVE N_G_DOF INDEX = N_G_DOF * (I_PT - 1) + J_PARM END IF END FUNCTION GET_PT_DOF_INDEX FUNCTION GET_INDEX_AT_PT (I_PT) RESULT (INDEX) ! * * * * * * * * * * * * * * * * * * * * * * * * * * ! DETERMINE DEGREES OF FREEDOM NUMBERS AT A NODE ! * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants ! for N_G_DOF, DOF_PT DOF_PT_SUM IMPLICIT NONE INTEGER, INTENT(IN) :: I_PT INTEGER :: INDEX (N_G_DOF) INTEGER :: EQ_PT, J ! implied loop ! DOF_PT = NUMBER OF DOF AT EACH NODE: (0:MAX_NP) ! DOF_PT_SUM = TOTAL NUMBER OF DOF, INCLUDING THIS NODE: (-1:MAX_NP) ! N_G_DOF = NUMBER OF PARAMETERS (DOF) PER NODE ! I_PT = SYSTEM NODE NUMBER ! INDEX = SYSTEM DOF NOS OF NODAL DOF IF ( DOF_VARY ) THEN ! VARIABLE DOF PER NODE ! EQ_IJ = DOF_PT_SUM (I-1) + J, J <= DOF_PT(I), I > 0 INDEX = 0 ! DEFAULT MAXIMUM ARRAY EQ_PT = DOF_PT (I_PT) ! DOF AT THIS NODE IF ( EQ_PT < 1 ) RETURN ! GEOMETRIC ONLY DO J = 1, EQ_PT INDEX (J) = DOF_PT_SUM (I_PT-1) + J END DO ! ON J ELSE ! ALL NODES HAVE N_G_DOF ! INDEX (J) = N_G_DOF * (I_PT - 1) + J INDEX = (/ (N_G_DOF * (I_PT - 1) + J, J = 1, N_G_DOF) /) END IF END FUNCTION GET_INDEX_AT_PT FUNCTION GET_ELEM_INDEX (LT_N, ELEM_NODES) RESULT(INDEX) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * ! DETERMINE DEGREES OF FREEDOM NUMBERS OF ELEMENT ! * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants ! for N_G_DOF, DOF_PT DOF_PT_SUM IMPLICIT NONE INTEGER, INTENT(IN) :: LT_N, ELEM_NODES (LT_N) INTEGER :: INDEX (LT_N * N_G_DOF) ! OUT INTEGER :: EQ_ELEM, EQ_PT, EQ_SYS, IG, K, SYS_K ! LOOPS ! INTEGER :: DOF_PT (0:MAX_NP), DOF_PT_SUM (-1:MAX_NP) ! DOF_PT = NUMBER OF DOF AT EACH NODE: (0:MAX_NP) ! DOF_PT_SUM = TOTAL NUMBER OF DOF, INCLUDING THIS NODE: (-1:MAX_NP) ! N_G_DOF = NUMBER OF GENERAL PARAMETERS (DOF) PER NODE ! LT_N = NUMBER OF NODES PER ELEMENT ! ELEM_NODES = NODAL INCIDENCES OF THE ELEMENT ! EQ_ELEM = LOCAL EQUATION NUMBER ! EQ_SYS = SYSTEM EQUATION NUMBER ! INDEX = SYSTEM DOF NUMBERS OF ELEMENT DOF NUMBERS IF ( DOF_VARY ) THEN ! VARIABLE DOF PER ELEMENT NODE ! EQ_IJ = DOF_PT_SUM (I-1) + J, J <= DOF_PT(I), I > 0 INDEX = 0 ! DEFAULT MAXIMUM ARRAY EQ_ELEM = 0 DO K = 1, LT_N ! LOOP OVER NODES OF ELEMENT SYS_K = ELEM_NODES (K) ! SYSTEM NODE NUMBER IF ( SYS_K < 1 ) CYCLE ! TO A REAL NODE EQ_PT = DOF_PT (SYS_K) ! DOF AT THIS NODE IF ( EQ_PT < 1 ) CYCLE ! TO A NON_GEOMETRIC NODE DO IG = 1, EQ_PT ! LOOP OVER GENERALIZED DOF EQ_ELEM = EQ_ELEM + 1 ! LOCAL EQ EQ_SYS = DOF_PT_SUM (SYS_K - 1) + IG ! SYSTEM EQ INDEX (EQ_ELEM) = EQ_SYS END DO ! OVER DOF END DO ! OVER LOCAL NODES ELSE ! ALL NODES HAVE N_G_DOF ! INDEX (N_G_DOF*(K-1)+IG) = N_G_DOF*(ELEM_NODES(K)-1) + IG DO K = 1, LT_N ! LOOP OVER NODES OF ELEMENT SYS_K = ELEM_NODES (K) ! SYSTEM NODE NUMBER DO IG = 1, N_G_DOF ! LOOP OVER GENERALIZED DOF EQ_ELEM = IG + N_G_DOF * (K - 1) ! LOCAL EQ EQ_SYS = IG + N_G_DOF * (SYS_K - 1) ! SYSTEM EQ IF ( SYS_K > 0 ) THEN ! VALID NODE INDEX (EQ_ELEM) = EQ_SYS ELSE ! ALLOW MISSING NODE INDEX (EQ_ELEM) = 0 END IF ! MISSING NODE END DO ! OVER DOF END DO ! OVER LOCAL NODES END IF ! VARIABLE DOF END FUNCTION GET_ELEM_INDEX !FUNCTION GET_ELEM_INDEX (LT_N, ELEM_NODES) RESULT(INDEX) !! * * * * * * * * * * * * * * * * * * * * * * * * * * !! DETERMINE DEGREES OF FREEDOM NUMBERS OF ELEMENT !! * * * * * * * * * * * * * * * * * * * * * * * * * * !Use System_Constants ! for N_G_DOF ! IMPLICIT NONE ! INTEGER, INTENT(IN) :: LT_N, ELEM_NODES (LT_N) ! INTEGER :: INDEX (LT_N * N_G_DOF) ! INTEGER :: IDOF, IELM, IG, K, CONST ! !! LT_N = NUMBER OF NODES PER ELEMENT !! N_G_DOF = NUMBER OF PARAMETERS (DOF) PER NODE !! ELEM_NODES = NODAL INCIDENCES OF THE ELEMENT !! INDEX = SYSTEM DOF NOS OF ELEMENT DOF ! ! DO K = 1, LT_N ! LOOP OVER NODES OF ELEMENT ! IDOF = - N_G_DOF ! IF ( ELEM_NODES (K) > 0) IDOF = IDOF + N_G_DOF * ELEM_NODES (K) ! CONST = N_G_DOF * (K - 1) !! LOOP OVER GENERALIZED DEGREES OF FREEDOM ! DO IG = 1, N_G_DOF ! IELM = CONST + IG !! INDEX (N_G_DOF * (K-1)+IG) = N_G_DOF * (ELEM_NODES(K)-1) + IG ! INDEX (IELM) = IDOF + IG ! END DO ! END DO !END FUNCTION GET_ELEM_INDEX FUNCTION GET_ELEM_NODES (L_ID, LT_N, NODES) RESULT (ELEM_NODES) ! * * * * * * * * * * * * * * * * * * * * * * * * ! EXTRACT NODES ASSOCIATED WITH ELEMENT L_ID ! * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants IMPLICIT NONE INTEGER, INTENT(IN) :: L_ID, LT_N INTEGER, INTENT(IN) :: NODES (L_S_TOT, NOD_PER_EL) INTEGER :: ELEM_NODES (LT_N) ! ELEM_NODES = THE NOD_PER_EL NODAL INCIDENCES OF THE ELEMENT ! L_ID = ELEMENT NUMBER ! LT_N = NUMBER OF NODES FOR THIS ELEMENT TYPE ! L_S_TOT = N_ELEMS + N_MIXED + N_F_SEG, TOTAL GENERALIZED ELEMENTS ! N_ELEMS = NUMBER OF ELEMENTS IN SYSTEM ! NODES = NODAL INCIDENCES OF ALL ELEMENTS ! NOD_PER_EL = NUMBER OF NODES PER ELEMENT ELEM_NODES (1:LT_N) = NODES (L_ID, 1:LT_N) if ( minval(elem_nodes) < 0 .or. & maxval(elem_nodes) > MAX_NP) & stop 'invalid elem_nodes, GET_ELEM_NODES' END FUNCTION GET_ELEM_NODES SUBROUTINE GET_ELEM_QUADRATURES ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! RECOVER QUADRATURE DATA BASED ON ELEMENT SHAPE ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use Precision_Module Use Elem_Type_Data ! for: ! LT_PARM, LT_QP, PT (LT_PARM, LT_QP), WT (LT_QP) IMPLICIT NONE INTEGER :: N_GP, NIP ! work, 1D rule Number of Integration Pts ! LT_SHAP = ELEMENT SHAPE CODE ! LT_QP = NUMBER OF QUADRATURE POINTS ! LT_PARM = DIMENSION OF SPACE ! PT = QUADRATURE COORDINATES ! WT = QUADRATURE WEIGHTS IF ( .NOT. TYPE_GAUS_ALLOC ) STOP 'ALLOC ERROR, GET_ELEM_QUADRATURES' N_GP = LT_QP SELECT CASE ( LT_SHAP ) CASE (0:1) ! LINE ELEMENT CALL GAUSS_1D (N_GP, N_GP, PT, WT) CASE (2) ! TRIANGULAR (UNIT COORDINATES) CALL DUNAVANT_UNIT_TRIANGLE_RULE (N_GP, PT, WT) CASE (3) ! QUADRILATERAL SELECT CASE ( LT_QP ) ! integer squared only CASE (0:1) ; NIP = 1 CASE (4) ; NIP = 2 CASE (9) ; NIP = 3 CASE (16) ; NIP = 4 CASE DEFAULT NIP = SQRT (FLOAT (N_GP) ) + 0.1 !9 ceiling END SELECT CALL GAUSS_2D (NIP, N_GP, PT, WT) CASE (4) ! HEXAHEDRA SELECT CASE ( LT_QP ) ! integer cubed only CASE (0:1) ; NIP = 1 CASE (8) ; NIP = 2 CASE (27) ; NIP = 3 CASE (64) ; NIP = 4 CASE DEFAULT NIP = (FLOAT (N_GP) ) ** (1. / 3.) + 0.1 !9 ceiling END SELECT CALL GAUSS_3D (NIP, N_GP, PT, WT) CASE (5) ! TETRAHEDRA (UNIT COORDINATES) CALL KEAST_UNIT_TET_RULE (N_GP, PT, WT) CASE (6) ! WEDGE STOP 'WEDGE NOT IN GET_ELEM_QUADRATURES' !b CASE (7) ! USER SUPPLIED QUADRATURE SUBROUTINE CASE DEFAULT ! UNSUPPORTED STOP 'INVALID ELEMENT SHAPE OPTION, GET_ELEM_QUADRATURES' END SELECT IF ( LT_QP == N_GP) RETURN WRITE (6, * ) 'LT_SHAP, LT_QP, NIP, N_GP', LT_SHAP, & LT_QP, NIP, N_GP, ' FATAL ERROR, GET_ELEM_QUADRATURES' STOP 'FATAL ERROR, GET_ELEM_QUADRATURES' END SUBROUTINE GET_ELEM_QUADRATURES SUBROUTINE GET_NUMB_QUADRATURES (PT_N, WT_N, N_GP) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! RECOVER QUADRATURE BASED ON NUMBER OF POINTS ! SEE: FUNCTION GET_NUMB_INERTIA_PTS ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use Precision_Module Use Elem_Type_Data ! for LT_PARM, LT_SHAP IMPLICIT NONE INTEGER, INTENT(IN) :: N_GP REAL(DP), INTENT(OUT) :: PT_N (LT_PARM, N_GP), WT_N (N_GP) INTEGER :: NIP ! LT_SHAP = ELEMENT SHAPE CODE ! N_GP = NUMBER OF QUADRATURE POINTS ! LT_PARM = DIMENSION OF SPACE ! PT_N = QUADRATURE COORDINATES ! WT_N = QUADRATURE WEIGHTS ! SELECT CASE ( LT_SHAP ) CASE (0:1) ! LINE ELEMENT CALL GAUSS_1D (N_GP, N_GP, PT_N, WT_N) CASE (2) ! TRIANGULAR (UNIT COORDINATES) CALL DUNAVANT_UNIT_TRIANGLE_RULE (N_GP, PT_N, WT_N) CASE (3) ! QUADRILATERAL NIP = SQRT (FLOAT (N_GP) ) + 0.1 CALL GAUSS_2D (NIP, N_GP, PT_N, WT_N) CASE (4) ! HEXAHEDRA NIP = (FLOAT (N_GP) ) ** (1. / 3.) + 0.1 CALL GAUSS_3D (NIP, N_GP, PT_N, WT_N) CASE (5) ! TETRAHEDRA (UNIT COORDINATES) CALL KEAST_UNIT_TET_RULE (N_GP, PT_N, WT_N) CASE (6) ! WEDGE STOP 'WEDGE NOT IN GET_NUMB_QUADRATURES' !b CASE (7) ! USER SUPPLIED QUADRATURE SUBROUTINE CASE DEFAULT ! UNSUPPORTED STOP 'INVALID ELEMENT SHAPE OPTION, GET_NUMB_QUADRATURES' END SELECT END SUBROUTINE GET_NUMB_QUADRATURES SUBROUTINE GET_ONE_PT_RULE (ONE_PT, ONE_WT) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! RECOVER ONE POINT QUADRATURE RULE ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use Precision_Module Use Elem_Type_Data ! for LT_PARM, LT_SHAP IMPLICIT NONE REAL(DP), INTENT(OUT) :: ONE_PT (LT_PARM), ONE_WT INTEGER, PARAMETER :: N_GP = 1 REAL(DP) :: PT_N (LT_PARM, N_GP), WT_N (N_GP) ! LT_SHAP = ELEMENT SHAPE CODE ! N_GP = NUMBER OF QUADRATURE POINTS ! LT_PARM = DIMENSION OF SPACE ! PT_N = QUADRATURE COORDINATES ! WT_N = QUADRATURE WEIGHTS ! SELECT CASE ( LT_SHAP ) CASE (0:1) ! LINE ELEMENT CALL GAUSS_1D (N_GP, N_GP, PT_N, WT_N) CASE (2) ! TRIANGULAR (UNIT COORDINATES) CALL DUNAVANT_UNIT_TRIANGLE_RULE (N_GP, PT_N, WT_N) CASE (3) ! QUADRILATERAL CALL GAUSS_2D (N_GP, N_GP, PT_N, WT_N) CASE (4) ! HEXAHEDRA CALL GAUSS_3D (N_GP, N_GP, PT_N, WT_N) CASE (5) ! TETRAHEDRA (UNIT COORDINATES) CALL KEAST_UNIT_TET_RULE (N_GP, PT_N, WT_N) CASE (6) ! WEDGE STOP 'WEDGE NOT IN GET_ONE_PT_RULE' !b CASE (7) ! USER SUPPLIED QUADRATURE SUBROUTINE CASE DEFAULT ! UNSUPPORTED STOP 'INVALID ELEMENT SHAPE OPTION, GET_ONE_PT_RULE' END SELECT ONE_PT (:) = PT_N (:, 1) ; ONE_WT = WT_N (1) ! copy END SUBROUTINE GET_ONE_PT_RULE !b FUNCTION GET_NUMB_INERTIA_PTS () RESULT (N_GP) !b ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * !b ! RECOVER NUMBER OF POINTS TO INTEGRATE INERTIA TENSOR !b ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * !b Use Precision_Module !b Use Elem_Type_Data !b IMPLICIT NONE !b INTEGER :: N_GP !b !b ! LT_SHAP = ELEMENT SHAPE CODE !b ! N_GP = NUMBER OF QUADRATURE POINTS !b ! LT_PARM = DIMENSION OF SPACE !b ! !b SELECT CASE ( LT_SHAP ) !b CASE (0:1) ; N_GP = 2 ! LINE ELEMENT !b CASE (2) ; N_GP = 3 ! TRIANGULAR ELEMENT !b CASE (3) ; N_GP = 4 ! QUADRILATERAL ELEMENT !b CASE (4) ; N_GP = 8 ! HEXAHEDRAL ELEMENT !b CASE (5) ; N_GP = 4 ! TETRAHEDRAL ELEMENT !b CASE (6) ; N_GP = 6 ! WEDGE ELEMENT !b CASE (7) ! USER SUPPLIED !b ! CALL USER_INERTIA_PT (N_GP, LT_PARM, PT, WT) !b !b STOP 'USER OPTION MISSING, GET_NUMB_INERTIA_PTS' !b CASE DEFAULT !b STOP 'INVALID OPTION, GET_NUMB_INERTIA_PTS' !b END SELECT !b END FUNCTION GET_NUMB_INERTIA_PTS SUBROUTINE GET_ELEM_TYPE_DATA (LT) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! GET DATA BASED ON ELEMENT TYPE (SEE CONTROL_INPUT ) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants Use Elem_Type_Data ! for LAST_LT, N_L_TYPE IMPLICIT NONE INTEGER, INTENT(IN) :: LT ! LAST_LT = LAST ACTIVE ELEMENT TYPE NUMBER ! LT = ELEMENT TYPE NUMBER ! LT_N = NUMBER OF NODES PER ELEMENT ! LT_QP = NUMBER OF QUADRATURE POINTS ! LT_GEOM = NUMBER OF ELEMENT GEOMERTY NODES ! LT_PARM = NUMBER OF PARAMETRIC SPACES FOR ELEMENT ! LT_FACES = NUMBER OF (FACING) NEIGHBORS ! LT_SONS = NUMBER OF SONS IF DIVIDED ! LT_SHAP = ELEMENT SHAPE FLAG NUMBER ! LT_USER = APPLICATION DEPENDENT OPTIONAL USER ITEM ! N_L_TYPE = NUMBER OF ELEMENT TYPES IF ( LT < 1 .OR. LT > N_L_TYPE ) STOP & 'ELEMENT TYPE IMPOSSIBLE, GET_ELEM_TYPE_DATA' ! GET CONTROLS FOR THIS TYPE LT_N = LT_DATA (1, LT) LT_QP = LT_DATA (2, LT) LT_GEOM = LT_DATA (3, LT) LT_PARM = LT_DATA (4, LT) LT_SHAP = LT_DATA (5, LT) LT_FACES = LT_DATA (6, LT) LT_SONS = LT_DATA (7, LT) LT_USER = LT_DATA (8, LT) LT_FREE = LT_N * N_G_DOF END SUBROUTINE GET_ELEM_TYPE_DATA SUBROUTINE LIST_ELEM_TYPE_DATA (LT) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! LIST DATA BASED ON ELEMENT TYPE (SEE CONTROL_INPUT ) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants Use Elem_Type_Data ! for LAST_LT, N_L_TYPE IMPLICIT NONE INTEGER, INTENT(IN) :: LT ! LAST_LT = LAST ACTIVE ELEMENT TYPE NUMBER ! LT = ELEMENT TYPE NUMBER ! LT_N = NUMBER OF NODES PER ELEMENT ! LT_QP = NUMBER OF QUADRATURE POINTS ! LT_GEOM = NUMBER OF ELEMENT GEOMERTY NODES ! LT_PARM = NUMBER OF PARAMETRIC SPACES FOR ELEMENT ! LT_FACES = NUMBER OF (FACING) NEIGHBORS ! LT_SONS = NUMBER OF SONS IF DIVIDED ! LT_SHAP = ELEMENT SHAPE FLAG NUMBER ! LT_USER = APPLICATION DEPENDENT OPTIONAL USER ITEM ! N_L_TYPE = NUMBER OF ELEMENT TYPES IF ( LT < 1 .OR. LT > N_L_TYPE ) THEN PRINT *, 'WARNING: ELEMENT TYPE ', & LT, ' IS IMPOSSIBLE, LIST_ELEM_TYPE_DATA' N_WARN = N_WARN + 1 ! INCREMENT WARNING END IF LT_FREE = LT_N * N_G_DOF ! LIST CONTROLS FOR THIS TYPE PRINT *, ' ' ; PRINT *, '*** LIST OF DATA FOR ELEMENT TYPE ', LT, ' ***' PRINT *, 'LT_N = ', LT_N PRINT *, 'LT_QP = ', LT_QP PRINT *, 'LT_GEOM = ', LT_GEOM PRINT *, 'LT_PARM = ', LT_PARM PRINT *, 'LT_SHAP = ', LT_SHAP PRINT *, 'LT_FACES = ', LT_FACES PRINT *, 'LT_SONS = ', LT_SONS PRINT *, 'LT_USER = ', LT_USER PRINT *, 'LT_FREE = ', LT_FREE END SUBROUTINE LIST_ELEM_TYPE_DATA SUBROUTINE SET_ELEM_TYPE_INFO (LT) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! SET DATA BASED ON ELEMENT TYPE (SEE CONTROL_INPUT ) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants Use Elem_Type_Data ! for: MAX_DAT, N_L_TYPE, & ! LT_DATA (MAX_DAT, N_L_TYPE), & ! LT_FREE, LT_GEOM, LT_N, LT_PARM, LT_QP, ELEM_NODES (LT_N), & ! COORD (LT_N, N_SPACE), GEOMETRY (LT_GEOM, N_SPACE), & ! C (LT_FREE), D (LT_FREE), S (LT_FREE, LT_FREE), & ! DLG (LT_PARM, LT_GEOM), DLG_QP (LT_PARM, LT_GEOM, LT_QP) & ! DLH (LT_PARM, LT_N), DLH_QP (LT_PARM, LT_N, LT_QP), & ! DLV (LT_PARM, LT_FREE), DLV_QP (LT_PARM, LT_FREE, LT_QP), & ! G (LT_GEOM), G_QP (LT_GEOM, LT_QP), & ! H (LT_N), H_QP (LT_N, LT_QP), & ! V (LT_FREE), V_QP (LT_FREE, LT_QP), & ! PT (LT_PARM, LT_QP), WT (LT_QP) IMPLICIT NONE INTEGER, INTENT(IN) :: LT ! LAST_LT = LAST ACTIVE ELEMENT TYPE NUMBER ! LT = ELEMENT TYPE NUMBER ! LT_DATA = ARRAY OF CONSTANTS FOR EACH ELEMENT TYPE ! 1, LT_N = NUMBER OF NODES PER ELEMENT ! 2, LT_QP = NUMBER OF QUADRATURE POINTS ! 3, LT_GEOM = NUMBER OF ELEMENT GEOMERTY NODES ! 4, LT_PARM = NUMBER OF PARAMETRIC SPACES FOR ELEMENT ! 5, LT_FACES = NUMBER OF (FACING) NEIGHBORS ! 6, LT_SONS = NUMBER OF SONS IF DIVIDED ! 7, LT_SHAP = ELEMENT SHAPE FLAG NUMBER ! 8, LT_USER = APPLICATION DEPENDENT OPTIONAL USER ITEM ! LT_FREE = LT_N * N_G_DOF ! N_L_TYPE = NUMBER OF ELEMENT TYPES IF ( LT < 1 .OR. LT > N_L_TYPE ) STOP & 'ELEMENT TYPE WRONG, SET_ELEM_TYPE_INFO' ! GET CONTROLS FOR THIS TYPE: ! LT_N, LT_QP, LT_GEOM, LT_PARM, ! LT_SHAP, LT_FACES, LT_SONS, LT_USER, LT_FREE CALL GET_ELEM_TYPE_DATA (LT) ! RESET LAST, ALLOCATE ARRAY FOR AN ELEMENT TYPE IF ( LT /= LAST_LT ) THEN LAST_LT = LT IF ( TYPE_APLY_ALLOC ) CALL DEALLOCATE_TYPE_APPLICATION CALL ALLOCATE_TYPE_APPLICATION IF ( TYPE_NTRP_ALLOC ) CALL DEALLOCATE_TYPE_INTERPOLATIONS CALL ALLOCATE_TYPE_INTERPOLATIONS !--> GET QUADRATURE RULE FOR ELEMENT TYPE AND SHAPE IF ( LT_QP > 0 ) THEN !b print *, 'SET_ELEM_TYPE_INFO LT, LT_QP ', LT, LT_QP !b CALL GET_ELEM_QUADRATURES CALL FILL_TYPE_INTERPOLATIONS END IF ! quadratures used END IF ! new element type END SUBROUTINE SET_ELEM_TYPE_INFO SUBROUTINE GET_REACTIONS (DD) ! * * * * * * * * * * * * * * * * * * * * * * * * * !--> COMPUTE REACTIONS IN FROM STORED SYSTEM ROWS ! * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants IMPLICIT NONE REAL(DP), INTENT(IN) :: DD (N_D_FRE) ! Automatic Arrays REAL(DP) :: TOTAL (N_G_DOF), SUM_POS (N_G_DOF), SUM_NEG (N_G_DOF) INTEGER :: IE, NODE, IG, J1, J2, J REAL(DP) :: R, CI, SIJ ! DD = COMPUTED SOLUTION VECTOR, SS*DD=CC ! MODE = STORAGE FLAG, 0=SYMMETRIC, OR 1=UNSYMMETRIC SKYLINE ! NOD_PER_EL = DOF NUMBER OF ESSENTIAL B.C. ! N_G_DOF = NUMBER OF DOF PER NODE ! N_REACT = SEQUENTIAL UNIT TO STORE REACTION DATA ! N_D_FRE = TOTAL NUMBER OF EQUATIONS ! TOTAL = SUM OF REACTIONS FOR N_G_DOF PARAMETERS IF ( N_REACT > 0 ) THEN REWIND N_REACT WRITE (6, 5) 5 FORMAT ( /, '*** REACTION RECOVERY ***', /, & ' NODE, PARAMETER, REACTION, EQUATION') TOTAL = 0.d0 ; SUM_POS = 0.d0 ; SUM_NEG = 0.d0 ELSE PRINT *,'WARNING: NO REACTION FILE, BUT GET_REACTIONS WAS CALLED' N_WARN = N_WARN + 1 RETURN END IF ! READ NODE, PARAMETER, RANGE OF NON-ZERO TERMS 10 READ (N_REACT, END = 30) NODE, IG, J1, J2 R = 0.d0 DO J = J1, J2 READ (N_REACT) SIJ R = R + SIJ * DD (J) END DO READ (N_REACT) CI R = R - CI IE = N_G_DOF * (NODE-1) + IG WRITE (6, 5010) NODE, DOF_NAME (IG), R, IE 5010 FORMAT ( I7, ', ', A12, 1X, 1PE11.4, I6 ) !b TOTAL (IG) = TOTAL (IG) + R IF ( R >= 0.d0 ) THEN !b SUM_POS (IG) = SUM_POS (IG) + R !b ELSE !b SUM_NEG (IG) = SUM_NEG (IG) + R !b END IF !b GO TO 10 !b replace with iostat in infinite loop !xxx 30 WRITE (6, 5020) 5020 FORMAT ('REACTION RESULTANTS', /, & 'PARAMETER, SUM POSITIVE NEGATIVE') TOTAL = SUM_POS + SUM_NEG !b DO J = 1, N_G_DOF WRITE (6, 5030) DOF_NAME (J), TOTAL (J), SUM_POS (J), SUM_NEG (J) 5030 FORMAT ( A12, 1PE11.4, 2X, 1PE11.4, 2X, 1PE11.4 ) END DO END SUBROUTINE GET_REACTIONS !SUBROUTINE INDXEL (NOD_PER_EL, N_EL_FRE, N_G_DOF, ELEM_NODES, INDEX) !! * * * * * * * * * * * * * * * * * * * * * * * * * * !! DETERMINE DEGREES OF FREEDOM NUMBERS OF ELEMENT !! * * * * * * * * * * * * * * * * * * * * * * * * * * !Use Precision_Module ! IMPLICIT NONE ! INTEGER, INTENT(IN) :: NOD_PER_EL, N_EL_FRE, N_G_DOF ! INTEGER, INTENT(IN) :: ELEM_NODES (NOD_PER_EL) ! INTEGER, INTENT(OUT) :: INDEX (N_EL_FRE) ! INTEGER :: CONST, IDOF, IG, IELM, K ! !! NOD_PER_EL = NUMBER OF NODES PER ELEMENT !! N_G_DOF = NUMBER OF PARAMETERS (DOF) PER NODE !! N_EL_FRE = NOD_PER_EL*N_G_DOF = NUMBER OF DOF PER ELEMENT !! ELEM_NODES = NODAL INCIDENCES OF THE ELEMENT !! INDEX = SYSTEM DOF NOS OF ELEMENT DOF ! ! !b IF ( DOF_VARY ) STOP 'INDXEL: INVALID FOR DOF_VARY' ! !! LOOP OVER NODES OF ELEMENT ! DO K = 1, NOD_PER_EL ! IDOF = - N_G_DOF ! IF ( ELEM_NODES (K) > 0) IDOF = IDOF + N_G_DOF * ELEM_NODES (K) ! CONST = N_G_DOF * (K - 1) !! LOOP OVER GENERALIZED DEGREES OF FREEDOM ! DO IG = 1, N_G_DOF ! IELM = CONST + IG !! INDEX (N_G_DOF * (K-1)+IG) = N_G_DOF * (ELEM_NODES(K)-1) + IG ! INDEX (IELM) = IDOF + IG ! END DO ! END DO !END SUBROUTINE INDXEL ! !SUBROUTINE INDXPT (I_PT, N_G_DOF, INDEX) !! * * * * * * * * * * * * * * * * * * * * * * * * * * !! DETERMINE DEGREES OF FREEDOM NUMBERS AT A NODE !! * * * * * * * * * * * * * * * * * * * * * * * * * * ! IMPLICIT NONE ! INTEGER, INTENT(IN) :: I_PT, N_G_DOF ! INTEGER, INTENT(OUT) :: INDEX (N_G_DOF) ! INTEGER :: J, CONST ! !! I_PT = SYSTEM NODE NUMBER !! N_G_DOF = NUMBER OF PARAMETERS (DOF) PER NODE !! INDEX = SYSTEM DOF NOS OF NODAL DOF ! ! !b IF ( DOF_VARY ) STOP 'INDXPT: INVALID FOR DOF_VARY' ! ! CONST = N_G_DOF * (I_PT - 1) ! DO J = 1, N_G_DOF !! INDEX (J) = N_G_DOF*(I_PT - 1) + J ! INDEX (J) = CONST + J ! END DO !END SUBROUTINE INDXPT SUBROUTINE ELEM_NODE_PROPERTIES (LT_N, PRT_L_PT, L_PT_PROP, & ELEM_NODES) ! * * * * * * * * * * * * * * * * * * * * * * * * * * ! EXTRACT FLOATING POINT PROPERTIES AT NODAL POINTS ! OF AN ELEMENT ! * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants Use Sys_Properties_Data IMPLICIT NONE INTEGER, INTENT(IN) :: LT_N INTEGER, INTENT(IN) :: ELEM_NODES (LT_N) REAL(DP), INTENT(OUT) :: PRT_L_PT (NOD_PER_EL, N_NP_FLO) INTEGER, INTENT(OUT) :: L_PT_PROP (NOD_PER_EL, N_NP_FIX) INTEGER :: I, IROW, J ! FLT_NP = FLOATING POINT PROP ARRAY OF SYSTEM NODES ! PRT_L_PT = FLOATING POINT PROP ARRAY OF ELEMENT NODES ! NP_FIX = INTEGER PROPERTY ARRAY OF SYSTEM NODES ! L_PT_PROP = INTEGER PROPERTY ARRAY OF ELEMENT NODES ! ELEM_NODES = ELEMENT INCIDENCES ARRAY OF THE ELEMENT ! MAX_NP = NUMBER OF SYSTEM NODES ! LT_N = NUMBER OF NODES PER ELEMENT ! N_NP_FIX = NUMBER OF INTEGER PROPERTIES PER NODE ! N_NP_FLO = NUMBER OF REAL PROPERTIES PER NODE ! N_HOMO = 1, IF PROPERTIES ARE SAME AT EACH NODE DO I = 1, LT_N IROW = ELEM_NODES (I) ! ALLOW FOR OMITTED NODES !9 use WHERE IF ( IROW > 0 ) THEN IF ( N_HOMO == 1) IROW = 1 IF ( N_NP_FLO > 0 ) THEN DO J = 1, N_NP_FLO PRT_L_PT (I, J) = FLT_NP (IROW, J) END DO END IF IF ( N_NP_FIX > 0 ) THEN DO J = 1, N_NP_FIX L_PT_PROP (I, J) = NP_FIX (IROW, J) END DO END IF END IF END DO END SUBROUTINE ELEM_NODE_PROPERTIES SUBROUTINE MATERIAL_PROP (NUM, PRT_MAT) !b not right ! * * * * * * * * * * * * * * * * * * * * * * * * * * ! EXTRACT REAL MATERIAL PROPERTIES OF MATERIAL ! NUM FROM MISCELLANEOUS REAL SYSTEM PROPERTIES ! * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants Use Sys_Properties_Data IMPLICIT NONE INTEGER, INTENT(IN) :: NUM REAL(DP), INTENT(OUT) :: PRT_MAT (MISC_FL) INTEGER :: I, ISTART ! NUM = MATERIAL NUMBER ! N_MAT = MAXIMUM ALLOWABLE MATERIAL NUMBER ! N_LP_FLO = NUMBER OF REAL ELEMENT PROPERTIES ! MISC_FL = NO OF MISC REAL SYSTEM PROPERTIES ! FLT_MISC = SYSTEM STORAGE FOR MISC REAL PROP ! PRT_MAT = REAL PROPERTY ARRAY FOR MATERIAL NUM !9 see if the following is still true ! PROPERTIES ARE STORED IN FLT_MISC IN ORDER OF MAT NO if ( num >= 0 ) stop 'MATERIAL_PROP invalid code' !b IF ( N_LP_FLO < 1 ) STOP 'BAD N_LP_FLO, MATERIAL_PROP' IF ( NUM > N_MAT ) STOP 'DIMENSIONS EXCEEDED, MATERIAL_PROP' ISTART = N_LP_FLO * (NUM - 1) DO I = 1, N_LP_FLO PRT_MAT (I) = FLT_MISC (ISTART + I) END DO END SUBROUTINE MATERIAL_PROP SUBROUTINE PENLTY (NP_FREE, COEF_C_EQ, C_Pen, S_Pen, WT) ! * * * * * * * * * * * * * * * * * * * * * * * * * * ! DEFINE CONSTRAINT PENALTY SQUARE AND COLUMN MATRICES ! * * * * * * * * * * * * * * * * * * * * * * * * * * Use Precision_Module INTEGER, INTENT(IN) :: NP_FREE REAL(DP), INTENT(IN) :: COEF_C_EQ (NP_FREE), WT REAL(DP), INTENT(OUT) :: C_Pen (NP_FREE), S_Pen (NP_FREE, NP_FREE) ! NP_FREE = NO DOF IN CONSTRAINT EQUATION ! COEF_C_EQ (I) = CONSTR EQ COEFFICIENT I+1 ! C_Pen = CONSTRAINT COLUMN MATRIX ! S_Pen = CONSTRAINT SQUARE MATRIX ! WT = PENALTY WEIGHT FACTOR ! INITIAL CALCULATIONS C_Pen (1) = 1.0 TEMP = COEF_C_EQ (NP_FREE) IF ( NP_FREE > 1 ) THEN DO I = 2, NP_FREE C_Pen (I) = COEF_C_EQ (I - 1) END DO END IF ! CALCULATE LEAST SQ CONSTRAINT FORMS DO I = 1, NP_FREE DO J = 1, NP_FREE S_Pen (J, I) = WT * C_Pen (I) * C_Pen (J) END DO END DO DO I = 1, NP_FREE C_Pen (I) = C_Pen (I) * TEMP * WT END DO END SUBROUTINE PENLTY SUBROUTINE POST_PROCESS_APPLICATION (NODES, DD, DD_OLD, ITER, X) !b SUBROUTINE POST_PROCESS_APPLICATION (NODES, DD, ITER, X) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! ELEMENT LEVEL POST-SOLUTION CALCULATIONS ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants Use Elem_Type_Data ! for PT (LT_PARM, LT_QP), WT (LT_QP), ! G (LT_GEOM, LT_QP), H (LT_N ), ! DLG (LT_PARM, LT_GEOM, LT_QP), ! DLH (LT_PARM, LT_N , LT_QP), ! S (LT_FREE, LT_FREE), C (LT_FREE), ! ELEM_NODES (LT_N) Use Sys_Properties_Data ! for FLT_MISC, MISC_FIX Use Select_Source Use Interface_Header IMPLICIT NONE ! Material option for future uses. Add to arguments, etc. ! ALWAYS USED REAL(DP), INTENT(INOUT) :: DD (N_D_FRE), DD_OLD (N_D_FRE) REAL(DP), INTENT(IN) :: X (MAX_NP, N_SPACE) INTEGER, INTENT(IN) :: NODES (L_S_TOT, NOD_PER_EL), ITER ! Automatic Arrays !b REAL(DP) :: EL_PROP (N_LP_FLO), PRT_MAT (MISC_FL) !b INTEGER :: L_PROP (N_LP_FIX), L_PT_PROP (N_NP_FIX) !b ?? !B REAL(DP) :: FIT_LT (LT_N, SCP_FIT) !b REAL(DP) :: PRT_MAT (MISC_FL) INTEGER :: L_PT_PROP (N_NP_FIX) !b ?? INTEGER :: IE, LT INTEGER, SAVE :: NOTE = 0 ! Allocatable Arrays REAL(DP), ALLOCATABLE :: PRT_L_PT (:, :) ! VARIABLES: ! COORD = COORDINATES OF ALL NODES ON ELEMENT ! D = NODAL PARAMETERS ASSOCIATED WITH AN ELEMENT ! DD = ARRAY OF SYSTEM DEGREES OF FREEDOM ! DGH = GLOBAL DERIVATIVES INTERPOLATION FUNCTIONS ! DLG = LOCAL DERIVATIVES GEOMETRIC INTERPOLATION ! DLH = LOCAL DERIVATIVES INTERPOLATION FUNCTIONS ! E = CONSTITUTIVE MATRIX ! EB = PRODUCT OF E*B ! FLT_MISC = SYSTEM STORAGE OF FLOATING PT MISC PROP ! FLUX_LT = FLUX AT ELEMENT NODES FROM A SCP: (SCP_FIT, LT_N) ! G = GEOMETRIC INTERPOLATION FUNCTIONS ! H = SOLUTION INTERPOLATION FUNCTIONS ! H_INTG = INTEGRAL OF INTERPOLATION FUNCTIONS ! INDEX = SYSTEM DOF NOS ASSOCIATED WITH ELEMENT ! ITER = CURRENT ITERATION NUMBER ! ELEM_NODES = THE NOD_PER_EL INCIDENCES OF THE ELEMENT ! L_PT_PROP = INTEGER PROPERTIES AT EACH ELEMENT NODE ! MAX_NP = TOTAL NUMBER OF NODES ! MISC_FIX = MISCELLANEOUS INTEGER SYSTEM PROPERTIES ! 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_EL_FRE = NUMBER OF DEGREES OF FREEDOM PER ELEMENT ! N_G_DOF = NUMBERS PARAMETERS PER NODE ! N_GEOM = NUMBER OF GEOMETRY NODES ! N_ITER = MAX NUMBER OF ITERATIONS ! N_MAT = NUMBER OF MATERIAL TYPES ! NODES = ELEMENT INCIDENCES OF ALL ELEMENTS ! N_PARM = DIMENSION OF PARAMETRIC SPACE ! N_QP = NUMBER OF QUADRATURE POINTS ! N_R_B = NUMBER OF ROWS IN B AND E MATRICES ! N_SPACE = DIMENSION OF SPACE ! N_FILE1 = UNIT FOR POST SOLUTION MATRICES STORAGE ! N_FILE2 = UNIT FOR POST SOLUTION MATRICES STORAGE ! N_FILE3-5 = OPTIONAL UNITS FOR USER (USED WHEN > 0) ! PRT_L_PT = REAL PROPERTIES AT ELEMENT NODES ! PRT_MAT = REAL ELEM PROPERTIES BASED ON MATERIAL NUMBER ! PT = QUADRATURE COORDINATES ! S = ELEMENT SQUARE MATRIX ! SCP_AVERAGES = AVERAGED FLUX ITEMS AT ALL NODES ! SCP_FIT = NUMBER IF TERMS BEING FIT, = N_R_B USUALLY ! STRAIN = STRAIN OR GRADIENT VECTOR ! STRAIN_0 = INITIAL STRAIN OR GRADIENT VECTOR ! STRESS = STRESS VECTOR ! WT = QUADRATURE WEIGHTS ! XYZ = SPACE COORDINATES AT A POINT ! N_FILE1 OR 2 MUST BE > 0 IF ( N_FILE1 == 0 .AND. N_FILE2 == 0 ) THEN STOP 'NO N_FILE1 OR NFILE2 IN POST_PROCESS_APPLICATION' ELSE IF ( N_FILE1 > 0 ) REWIND (N_FILE1) IF ( N_FILE2 > 0 ) REWIND (N_FILE2) END IF IF ( U_FLUX > 0 ) THEN REWIND (U_FLUX ) ELSE PRINT *, 'NOTE: UNIT U_FLUX NOT AVAILABLE FOR APPLICATION' END IF WRITE (N_PRT, "(/,'** BEGIN ELEMENT APPLICATION POST PROCESSING **')") LT = 1 ; LAST_LT = 0 ! INITIALIZE ELEMENT TYPES !b STRAIN = 0.d0 ; STRAIN_0 = 0.d0 ; STRESS = 0.d0 ! INITRIALIZE !b !--> LOOP OVER ELEMENTS DO IE = 1, L_S_TOT ! for elements and any boundary segments CALL SET_TO_MIXED_OR_SEG_OR_STD (IE) ! For IS_STD, IS_MIX, IS_SEG IF ( IE == (N_ELEMS + 1) ) WRITE (N_PRT, & "('** BEGIN APPLICATION BOUNDARY SEGMENTS POST PROCESSING **')") !--> 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) !--> Allocate Arrays for this type IF ( ALLOCATED (PRT_L_PT) ) DEALLOCATE ( PRT_L_PT ) !ba ALLOCATE ( PRT_L_PT ( LT_N, N_NP_FLO ) ) !b need the same for integer prop !b IF ( I_BUG > 0 ) PRINT *, ALLOCATED ( PRT_L_PT ) !--> GET QUADRATURE RULE FOR ELEMENT TYPE AND SHAPE IF ( LT_QP > 0 ) CALL GET_ELEM_QUADRATURES 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 !--> EXTRACT NODAL COORDINATES CALL ELEM_COORD (LT_N, N_SPACE, X, COORD, ELEM_NODES) !--> CALCULATE DEGREE OF FREEDOM NUMBERS INDEX = GET_ELEM_INDEX (LT_N, ELEM_NODES) !--> EXTRACT NODAL PARAMETERS OF THE ELEMENT D = GET_ELEM_DOF (DD) !--> EXTRACT NODAL POINT PROPERTIES (IF ANY) IF ( N_NP_FLO > 0) CALL ELEM_NODE_PROPERTIES (LT_N, & PRT_L_PT, L_PT_PROP, ELEM_NODES) !b !--> EXTRACT ELEMENT PROPERTIES (IF ANY) !b IF ( LP_TEST > 0) CALL ELEM_PROPERTIES (IE, L_PROP, EL_PROP) !B !--> EXTRACT ELEMENT NODAL FLUXES AT ELEM_NODES (IF ANY) !b !B IF ( SCP_RECORDS_SAVED ) THEN !b !B FIT_LT = GET_SCP_ELEM_FIT (SCP_AVERAGES) !b !B ELSE ; FIT_LT = 0.d0 ; END IF !b !--> EXTRACT MATERIAL PROPERTIES (IF ANY) IF ( N_MAT > 0) CALL MATERIAL_PROP (N_MAT, PRT_MAT) ! !--> PERFORM PROBLEM DEPENDENT CALCULATIONS AND OUTPUT ! !b IF ( POST_EL .AND. IS_ELEMENT ) THEN !b IF ( IE <= N_ELEMS ) THEN !b IF ( DEBUG_POST_EL ) PRINT *,'Debug IS_STD ', IS_STD IF ( POST_EL .AND. IS_ELEMENT ) CALL & SELECT_POST_PROCESS_ELEM (PRT_L_PT, PRT_MAT, & L_PT_PROP, ITER, IE, & DD, DD_OLD) !b !b L_PT_PROP, ITER, IE) !b ELSE IF ( POST_MIXED .AND. IS_MIXED_BC ) THEN !b IF ( DEBUG_MIX_SQ ) PRINT *,'Debug IS_MIX ', IS_MIX CALL SELECT_POST_PROCESS_MIXED (PRT_L_PT, PRT_MAT, & !b L_PT_PROP, ITER, IS_MIX) !b ELSE !b IF ( NOTE == 0 ) THEN ; NOTE = 1; PRINT *,'WARNING, SELECT_POST_PROCESS UNKNOWN OPTION' !b PRINT *, THIS_EL, 'NOT A STANDARD OR MIXED SEGMENT ELEMENT' PRINT *,'IS_STD, IS_MIX, IS_SEG', IS_STD, IS_MIX, IS_SEG PRINT *,'POST_EL, IS_ELEMENT', POST_EL, IS_ELEMENT PRINT *,'POST_MIXED, IS_MIXED_BC', POST_MIXED, IS_MIXED_BC N_WARN = N_WARN + 1 !b END IF ! Warn only once END IF ! Element class !b IF ( IE == N_ELEMS ) WRITE (N_PRT, & "('ELEMENT POST PROCESSING COMPLETE.')") IF ( IE == L_S_TOT .AND. (N_MIXED > 0 .OR. N_F_SEG > 0) ) & WRITE (N_PRT, "('POST PROCESSING COMPLETE.', /)") END DO ! over all elements and boundary segments !--> POST-PROCESSING COMPLETE !b CALL DEALLOCATE_TYPE_INTERPOLATIONS IF ( ALLOCATED (PRT_L_PT) ) DEALLOCATE ( PRT_L_PT ) END SUBROUTINE POST_PROCESS_APPLICATION SUBROUTINE POST_PROCESS_GRADS (NODES, DD, ITER) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! SAVE ELEMENT GRADIENTS AS SCP INPUT RECORDS ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants ! for L_S_TOT, N_D_FRE, NOD_PER_EL ! N_L_TYPE, THIS_EL, U_FLUX Use Elem_Type_Data ! for PT (LT_PARM, LT_QP), WT (LT_QP), ! G (LT_GEOM, LT_QP), H (LT_N ), ! DLG (LT_PARM, LT_GEOM, LT_QP), ! DLH (LT_PARM, LT_N , LT_QP), ! S (LT_FREE, LT_FREE), C (LT_FREE), ! ELEM_NODES (LT_N), D (LT_FREE) !b Use Sys_Properties_Data ! for FLT_MISC, MISC_FIX Use Interface_Header IMPLICIT NONE REAL(DP), INTENT(IN) :: DD (N_D_FRE) INTEGER, INTENT(IN) :: NODES (L_S_TOT, NOD_PER_EL), ITER INTEGER :: IE, LT ! VARIABLES: ! D = NODAL PARAMETERS ASSOCIATED WITH AN ELEMENT ! DD = ARRAY OF SYSTEM DEGREES OF FREEDOM ! DGH = GLOBAL DERIVATIVES INTERPOLATION FUNCTIONS ! E = CONSTITUTIVE MATRIX ! INDEX = SYSTEM DOF NOS ASSOCIATED WITH ELEMENT ! ITER = CURRENT ITERATION NUMBER ! 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_L_TYPE = NUMBER OF DIFFERENT ELEMENT TYPES USED ! N_ELEMS = NUMBER OF ELEMENTS IN SYSTEM ! NODES = ELEMENT INCIDENCES OF ALL ELEMENTS ! N_QP = NUMBER OF QUADRATURE POINTS ! N_R_B = NUMBER OF ROWS IN B AND E MATRICES ! U_FLUX = BINARY UNIT TO STORE GRADIENTS OR FLUXES ! U_FLUX MUST BE > 0 IF ( U_FLUX <= 0 ) THEN PRINT *, 'NO U_FLUX IN POST_PROCESS_GRADS' ELSE REWIND (U_FLUX ) END IF WRITE (N_PRT, "(/, 'BEGIN ELEMENT FLUX SAVE, ITER =', I4)") ITER LT = 1 ! INITIALIZE ELEMENT TYPES !b STRAIN = 0.d0 ; STRAIN_0 = 0.d0 ; STRESS = 0.d0 ! INITRIALIZE !b !b LAST_LT = 0 !b !b call inquire_about_file ("GRAD.BIN") !b !b call inquire_about_file ("N_FILE1.BIN") !b !--> LOOP OVER ELEMENTS DO IE = 1, N_ELEMS ! for elements and any boundary segments CALL SET_THIS_ELEMENT_NUMBER (IE) !--> 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) !B set or get ??? !b END IF ! a new element type ! INITIALIZE AND CHECK IF ACTIVE !b IF ( DELETED (IE) ) CYCLE ! to skip deleted element !XXX !--> 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) !--> EXTRACT NODAL PARAMETERS OF THE ELEMENT D = GET_ELEM_DOF (DD) !--> RECOVER FLUXES, LIST, SAVE FOR SCP !b print *, 'IE, LT, LT_QP, U_FLUX ', IE, LT, LT_QP, U_FLUX !b call LIST_SYSTEM_UNIT_NUMBERS !b IF ( USE_EXACT_FLUX ) THEN CALL LIST_ELEM_AND_EXACT_FLUXES (U_FLUX, IE) !b ELSE CALL LIST_ELEM_FLUXES (U_FLUX, IE) !b END IF IF ( IE == N_ELEMS ) WRITE (N_PRT, & "('ELEMENT FLUX SAVE COMPLETE (TO UNIT U_FLUX).')") END DO ! over all elements END SUBROUTINE POST_PROCESS_GRADS SUBROUTINE PT_COORD (I_PT, P_SPACE, X, COORD) ! * * * * * * * * * * * * * * * * * * * * * * * * * ! EXTRACT COORDINATES OF POINT NUMBER I_PT ! * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants INTEGER, INTENT(IN) :: I_PT, P_SPACE REAL(DP), INTENT(IN) :: X (MAX_NP, P_SPACE) REAL(DP), INTENT(OUT) :: COORD (1, P_SPACE) ! X = SPATIAL COORDINATES OF ALL SYSTEM NODES ! COORD = SPATIAL COORDINATES OF THE NODE ! MAX_NP = TOTAL NUMBER OF NODES IN SYSTEM ! NOD_PER_EL = NUMBER OF NODES PER ELEMENT ! P_SPACE = DIMENSION OF THE SPACE COORD (1, 1:P_SPACE) = X (I_PT, 1:P_SPACE) END SUBROUTINE PT_COORD FUNCTION GET_PT_COORD (I_PT, X) RESULT (COORD_PT) ! * * * * * * * * * * * * * * * * * * * * * * * * * ! EXTRACT COORDINATES OF POINT NUMBER I_PT ! * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants INTEGER, INTENT(IN) :: I_PT REAL(DP), INTENT(IN) :: X (MAX_NP, N_SPACE) REAL(DP) :: COORD_PT (N_SPACE) ! X = SPATIAL COORDINATES OF ALL SYSTEM NODES ! COORD_PT = SPATIAL COORDINATES OF THE NODE ! MAX_NP = TOTAL NUMBER OF NODES IN SYSTEM ! NOD_PER_EL = NUMBER OF NODES PER ELEMENT ! N_SPACE = DIMENSION OF THE SPACE COORD_PT (1:N_SPACE) = X (I_PT, 1:N_SPACE) END FUNCTION GET_PT_COORD SUBROUTINE START_ITERATION (X, DD) ! * * * * * * * * * * * * * * * * * * * * * * * * * * ! INITIALIZE SYSTEM DOF FOR ITERATIVE SOLUTION ! * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants Use Interface_Header IMPLICIT NONE REAL(DP), INTENT(IN) :: X (MAX_NP, N_SPACE) REAL(DP), INTENT(OUT) :: DD (N_D_FRE) REAL(DP) :: START_DOF_VALUE INTEGER :: I, INDX, J, IPRINT = 1 ! Automatic Arrays REAL(DP) :: XYZ_P (N_SPACE) INTEGER :: INDEX (N_G_DOF) ! ! MAX_NP = NUMBER OF NODES IN SYSTEM ! N_G_DOF = NUMBER OF PARAMETERS (DOF) PER NODE ! N_SPACE = DIMENSION OF SPACE ! N_D_FRE = TOTAL NUMBER OF SYSTEM DOF ! INDEX = SYSTEM DOF NUMBERS FOR DOF AT NODE ! X = COORDINATES OF SYSTEM NODES ! XYZ_P = SPATIAL COORDINATE ARRAY OF A NODE ! DD = SYSTEM ARRAY OF DEGREES OF FREEDOM ! IPRINT > 0, PRINT THE STARTING VALUES IF ( ECHO_START ) WRITE (N_PRT, 5) N_G_DOF ; 5 FORMAT ( /, & '** STARTING VALUES FOR ITERATIVE SOLUTION **',/, & 'NODE ', I2,' PARAMETER & VALUE PAIRS') DO I = 1, MAX_NP ! FIND PT COORDS AND DOF NOS INDEX = GET_INDEX_AT_PT (I) XYZ_P = GET_PT_COORD (I, X) IF ( ECHO_START ) WRITE (N_PRT, '(I5)', ADVANCE="NO") I DO J = 1, N_G_DOF INDX = INDEX (J) DD (INDX) = START_DOF_VALUE (J, XYZ_P) ! START IS A FUNCTION TO DEFINE INITIAL VALUES ! OF THE SYSTEM DEGREES OF FREEDOM IF ( ECHO_START ) WRITE (N_PRT, 510, ADVANCE="NO") J, DD (INDX) 510 FORMAT ( I4, 1PE13.5 ) END DO IF ( ECHO_START ) WRITE (N_PRT, '(" ")') END DO END SUBROUTINE START_ITERATION SUBROUTINE SET_INITIAL_CONDITION (X, DD) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * ! INITIALIZE SYSTEM DOF FOR TIME DEPENDENT INTEGRATION ! Keyword initial_function should be true ! * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants Use Interface_Header IMPLICIT NONE REAL(DP), INTENT(IN) :: X (MAX_NP, N_SPACE) REAL(DP), INTENT(OUT) :: DD (N_D_FRE) REAL(DP) :: START_DOF_VALUE INTEGER :: I, INDX, J, IPRINT = 1 ! Automatic Arrays REAL(DP) :: XYZ_P (N_SPACE) INTEGER :: INDEX (N_G_DOF) ! ! MAX_NP = NUMBER OF NODES IN SYSTEM ! N_G_DOF = NUMBER OF PARAMETERS (DOF) PER NODE ! N_SPACE = DIMENSION OF SPACE ! N_D_FRE = TOTAL NUMBER OF SYSTEM DOF ! INDEX = SYSTEM DOF NUMBERS FOR DOF AT NODE ! X = COORDINATES OF SYSTEM NODES ! XYZ_P = SPATIAL COORDINATE ARRAY OF A NODE ! DD = SYSTEM ARRAY OF DEGREES OF FREEDOM ! IPRINT > 0, PRINT THE STARTING VALUES IF ( IPRINT > 0 ) WRITE (N_PRT, 5) N_G_DOF 5 FORMAT ( /, '** SETTING DOF INITIAL CONDITIONS **',/, & 'NODE ', I2,' PARAMETER & VALUE PAIRS') DO I = 1, MAX_NP ! FIND PT COORDS AND DOF NOS INDEX = GET_INDEX_AT_PT (I) XYZ_P = GET_PT_COORD (I, X) IF ( IPRINT > 0 ) WRITE (N_PRT, '(I5)', ADVANCE="NO") I DO J = 1, N_G_DOF INDX = INDEX (J) DD (INDX) = START_DOF_VALUE (J, XYZ_P) !b DD (INDX) = SET_INITIAL_VALUE (J, XYZ_P) !b ! START IS A FUNCTION TO DEFINE INITIAL VALUES ! OF THE SYSTEM DEGREES OF FREEDOM IF ( IPRINT > 0 ) WRITE (N_PRT, 510, ADVANCE="NO") J, DD (INDX) 510 FORMAT ( I4, 1PE13.5 ) END DO IF ( IPRINT > 0 ) WRITE (N_PRT, '(" ")') END DO END SUBROUTINE SET_INITIAL_CONDITION SUBROUTINE STORE_COLUMN (LT_FREE, INDEX, C, CC) ! * * * * * * * * * * * * * * * * * * * * * * * * * * ! STORE ELEMENT COLUMN MATRIX IN SYSTEM COLUMN MATRIX ! * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants IMPLICIT NONE INTEGER, INTENT(IN) :: LT_FREE, INDEX (LT_FREE) REAL(DP), INTENT(IN) :: C (LT_FREE) REAL(DP), INTENT(INOUT) :: CC (N_D_FRE) INTEGER :: I, J ! INDEX = SYSTEM DOF NOS OF THE ELEMENT DOF ! C = ELEMENT COLUMN MATRIX ! CC = SYSTEM COLUMN MATRIX ! N_D_FRE = NO DEGREES OF FREEDOM IN THE SYSTEM ! LT_FREE = NUMBER OF DEGREES OF FREEDOM PER ELEMENT DO I = 1, LT_FREE J = INDEX (I) IF ( J > 0) CC (J) = CC (J) + C (I) END DO END SUBROUTINE STORE_COLUMN SUBROUTINE SUM_RHS (CC) ! * * * * * * * * * * * * * * * * * * * * * * * * * * ! SUM INPUT VALUES IN FORCING VECTOR, CC ! * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants IMPLICIT NONE REAL(DP), INTENT(IN) :: CC (N_D_FRE) ! Automatic Arrays REAL(DP) :: TOTAL (N_G_DOF), SUM_POS (N_G_DOF), SUM_NEG (N_G_DOF) !b INTEGER :: I, J, IJ_EQ ! N_D_FRE = TOTAL NUMBER OF SYSTEM DEGREES OF FREEDOM ! N_G_DOF = NUMBER OF PARAMETERS PER NODE ! CC = SYSTEM EQUATIONS COLUMN MATRIX ! MAX_NP = TOTAL NUMBER OF NODES IN THE SYSTEM TOTAL = 0.d0 ; SUM_POS = 0.d0 ; SUM_NEG = 0.d0 !b WRITE (6, 5030) 5030 FORMAT (/, '*** INPUT SOURCE RESULTANTS ***',/, & 'ITEM SUM POSITIVE NEGATIVE') DO I = 1, MAX_NP DO J = 1, N_G_DOF IJ_EQ = N_G_DOF * (I - 1) + J !b TOTAL (J) = TOTAL (J) + CC (IJ_EQ) IF ( CC (IJ_EQ) >= 0.d0 ) THEN !b SUM_POS (J) = SUM_POS (J) + CC (IJ_EQ) !b ELSE !b SUM_NEG (J) = SUM_NEG (J) + CC (IJ_EQ) !b END IF !b END DO END DO TOTAL = SUM_POS + SUM_NEG DO J = 1, N_G_DOF !b WRITE (6, '( I5, 2X, 1PE12.4 )') J, TOTAL (J) WRITE (6, '( I5, 3(2X, 1PE12.4) )') & J, TOTAL (J), SUM_POS (J), SUM_NEG (J) END DO END SUBROUTINE SUM_RHS SUBROUTINE CONVECTION_MATRICES (NF, NOD_PER_EL, C, S, & COORD, N_SPACE, IOPT, FLUX, NQ, PT, WT) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! FACE FLUX MATRICES ON PLANAR FACES (OR ELEMENTS) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use Precision_Module IMPLICIT NONE INTEGER, PARAMETER :: N_PARM = 2 ! ALWAYS USED INTEGER, INTENT(IN) :: NF, NOD_PER_EL, N_SPACE, NQ INTEGER, INTENT(OUT) :: IOPT REAL(DP), INTENT(IN) :: COORD (NOD_PER_EL, N_SPACE), FLUX (NF) REAL(DP), INTENT(IN) :: PT (N_PARM, NQ), WT (NQ) REAL(DP), INTENT(OUT) :: C (NOD_PER_EL) REAL(DP), INTENT(OUT) :: S (NOD_PER_EL, NOD_PER_EL) INTEGER :: I, J, IQ REAL(DP) :: F1, F2, F12, DET, DET_WT ! OPTIONAL FOR NUMERICAL INTEGRATION REAL(DP) :: H (NOD_PER_EL ), & DLH (N_SPACE, NOD_PER_EL ), & AJ (N_SPACE, N_SPACE), AJ_INV (N_SPACE, N_SPACE) ! AJ = JACOBIAN ! AJ_INV = JACOBIAN INVERSE ! C = BOUNDARY FLUX COLUMN MATRIX CONTRIBUTIONS ! COORD = SPATIAL COORDINATES OF SEGMENT NODES ! DLG = LOCAL DERIVATIVES GEOMETRIC INTERPOLATION ! DLH = LOCAL DERIVATIVES INTERPOLATION FUNCTIONS ! FLUX = SPECIFIED BOUNDARY FLUX COMPONENTS ! H = SOLUTION INTERPOLATION FUNCTIONS ! IOPT = PROBLEM MATRIX REQUIREMENTS, MUST BE SET. ! = 0, NO ACTION RETURN ! = 1, CALCULATE C ONLY ! = 2, CALCULATE S ONLY ! = 3, CALCULATE BOTH C AND S ! NOD_PER_EL = NO. OF NODES ON AN ELEMENT BOUNDARY SEGMENT ! NF = NUMBER OF FLUX COMPONENTS PER NODE ! N_D_FLUX = N*NF = MAXIMUM NUMBER OF FLUX CONTRIBUTIONS ! N_G_DOF = NUMBER OF PARAMETERS PER NODE POINT = ! ! N_PARM = PARAMETRIC GEOMETRY SPACE, = N_SPACE USUALLY ! NQ = NUMBER OF QUADRATURE POINTS ! N_SPACE = DIMENSION OF SOLUTION SPACE ! PT = QUADRATURE COORDINATES ! S = BOUNDARY FLUX SQUARE MATRIX ! WT = QUADRATURE WEIGHTS IOPT = 3 ! S & C computed F1 = FLUX (1) ! convection coefficient F2 = FLUX (2) ! reference temperature F12 = F1 * F2 ! COMPUTE MATRICES BY NUMERICAL INTEGRATION (ISOPARAMETRIC) C = 0.d0 ; S = 0.d0 ! ZERO WORKSPACES DO IQ = 1, NQ CALL SCALAR_SHAPES (PT (:, IQ), H) CALL SCALAR_DERIVS (PT (:, IQ), DLH) AJ = MATMUL ( DLH, COORD ) CALL INVERT_JACOBIAN (AJ, AJ_INV, DET, N_SPACE) DET_WT = ABS (DET) * WT (IQ) ! PLANAR FACE (THUS ELEMENT) CONVECTION DO I = 1, NOD_PER_EL C (I) = C (I) + H (I) * F12 * DET_WT DO J = 1, NOD_PER_EL S (I, J) = S (I, J) + H (I) * H (J) * DET_WT * F1 END DO END DO END DO END SUBROUTINE CONVECTION_MATRICES FUNCTION INVERSION_OF (A_Given) RESULT (A_inv) ! * * * * * * * * * * * * * * * * * * * * * * * * * * ! INVERSION OF NONSYMMETRIC MATRIX A_Given (N_size, N_size) ! * * * * * * * * * * * * * * * * * * * * * * * * * * Use Precision_Module IMPLICIT NONE REAL(DP), INTENT(IN) :: A_Given (:, :) ! Automatic Arrays REAL(DP) :: A_inv ( Size (A_Given, 1), Size (A_Given, 1) ) REAL(DP) :: B ( Size (A_Given, 1) ), C ( Size (A_Given, 1) ), D INTEGER :: I, J, K, M, N_size ! A IS DESTROYED AND REPLACED BY ITS INVERSE ! B,C = ARE WORKING SPACE VECTORS ! N_size = SIZE OF GIVEN MATRIX N_size = Size (A_Given, 1) A_inv = A_Given ! Copy into inverse IF ( A_inv (1, 1) == 0.0 ) STOP 'ZERO PIVOT IN INVERSION_OF' A_inv (1, 1) = 1.d0 / A_inv (1, 1) ! NOTE: several loops could be replaced with matmul DO M = 1, (N_size - 1) K = M + 1 DO I = 1, M B (I) = 0.d0 DO J = 1, M B (I) = B (I) + A_inv (I, J) * A_inv (J, K) END DO END DO D = 0.d0 DO I = 1, M D = D + A_inv (K, I) * B (I) END DO D = A_inv (K, K) - D IF ( D == 0.d0 ) STOP 'ZERO PIVOT IN INVERSION_OF' A_inv (K, K) = 1.d0 / D DO I = 1, M A_inv (I, K) = - B (I) * A_inv (K, K) END DO DO J = 1, M C (J) = 0.d0 DO I = 1, M C (J) = C (J) + A_inv (K, I) * A_inv (I, J) END DO END DO DO J = 1, M A_inv (K, J) = - C (J) * A_inv (K, K) END DO DO I = 1, M DO J = 1, M A_inv (I, J) = A_inv (I, J) - B (I) * A_inv (K, J) END DO END DO END DO END FUNCTION INVERSION_OF !SUBROUTINE REACTEL_ (IE, NOD_PER_EL, N_EL_FRE, N_G_DOF, N_FILE, & ! ELEM_NODES, IOPT, D) !! * * * * * * * * * * * * * * * * * * * * * * * * * * * !! GET REACTIONS (FLUXES) AT AN ELEMENTS NODES !! * * * * * * * * * * * * * * * * * * * * * * * * * * * !Use Precision_Module ! IMPLICIT NONE ! INTEGER, INTENT(IN) :: IE, NOD_PER_EL, N_EL_FRE, N_G_DOF ! INTEGER, INTENT(IN) :: N_FILE, IOPT ! INTEGER, INTENT(IN) :: ELEM_NODES (NOD_PER_EL) ! REAL(DP), INTENT(IN) :: D (N_EL_FRE) ! ! INTEGER :: IG, IN, IROW, L ! REAL(DP) :: S (N_EL_FRE, N_EL_FRE), C (N_EL_FRE) ! REAL(DP) :: SUM_PT (N_G_DOF), REACT, ROW ! REAL(DP) :: SUM_EL (N_G_DOF, NOD_PER_EL) ! !! C = ELEMENT COLUMN MATRIX !! D = KNOWN SOLUTION CAUSING THE REACTIONS !! IE = CURRENT ELEMENT NUMBER !! IOPT = WRITE OPTION CODE 1-ONCE, 2-TWICE !! ELEM_NODES = ELEMENT TOPOLOGY LIST !! NOD_PER_EL = NUMBER OF NODES PER ELEMENT !! N_EL_FRE = 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_EL = ELEMENT WORKSPACE RESULT !! SUM_PT = NODAL WORKSPACE ESULT ! ! IF ( N_FILE < 1) STOP 'INVALID UNIT IN REACTEL' !! GIVE REACTIONS, R = S*D - C ! SUM_PT = 0. ; SUM_EL = 0. !! READ STIFFNESS, AND SOURCE FOR ELEM REACTIONS !! SINGLE READ ! IF ( IOPT == 1 ) THEN ! READ (N_FILE) S, C !! DOUBLE READ ! ELSEIF (IOPT.EQ.2) THEN ! READ (N_FILE) S ; READ (N_FILE) C ! END IF ! DO IN = 1, NOD_PER_EL ! IF (IN == 1) THEN ! WRITE (6, * ) ' ELEMENT', IE, ' REACTIONS' ! WRITE (6, * ) ' NODE IDOF REACTION SOURCE' ! END IF ! DO IG = 1, N_G_DOF ! IROW = N_G_DOF * (IN - 1) + IG ! SUM_EL (IG, 1) = SUM_EL (IG, 1) + C (IROW) ! ROW = 0.D0 ! DO L = 1, N_EL_FRE ! ROW = ROW + S (IROW, L) * D (L) !9 use matmul ! END DO ! REACT = ROW - C (IROW) ! SUM_PT (IG) = SUM_PT (IG) + REACT ! PRINT "( I5, I5, 1PE15.5, 1PE15.5 )", & ! ELEM_NODES (IN), IG, REACT, C (IROW) ! END DO ! END DO ! DO IG = 1, N_G_DOF ! PRINT "(' SUM:', I5, 1PE15.5, 1PE15.5 )", & ! IG, SUM_PT (IG), SUM_EL (IG, 1) ! END DO !END SUBROUTINE REACTEL_ SUBROUTINE CAL_PRT (H, PRT_L_PT, VALUES) ! * * * * * * * * * * * * * * * * * * * * * * * * * * ! CALCULATE N_NP_FLO PROPERTIES AT A LOCAL PT USING ! ELEMENT'S NODAL PROPERTIES, PRT_L_PT, AND THE NOD_PER_EL ! INTERPOLATION FUNCTIONS, H, AT THE POINT ! * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants IMPLICIT NONE REAL(DP), INTENT(IN) :: H (NOD_PER_EL) REAL(DP), INTENT(IN) :: PRT_L_PT (NOD_PER_EL, N_NP_FLO) REAL(DP), INTENT(OUT) :: VALUES (N_NP_FLO) INTEGER :: I, J REAL(DP) :: SUM ! NOD_PER_EL = NUMBER OF NODES PER ELEMENT ! N_NP_FLO = NO. OF FLOATING POINT NODAL PROPERTIES ! H = C^0 INTERPOLATION FUNCTIONS FOR AN ELEMENT ! PRT_L_PT = FLOATING PT PROPS OF ELEMENT'S NODES ! VALUES = LOCAL VALUES OF PROPERTIES IF ( N_NP_FLO < 1) STOP 'N_NP_FLO = 0, CAL_PRT' DO I = 1, N_NP_FLO SUM = 0.0 DO J = 1, NOD_PER_EL SUM = SUM + H (J) * PRT_L_PT (J, I) END DO VALUES (I) = SUM END DO !9 use VALUES (I) = SUM ( H(:)*PRT_L_PT (:, I) ) END SUBROUTINE CAL_PRT ! SUBROUTINE GET_MATERIAL_PROPS (NUM, N_PR_MAT, FLT_MISC, PRT_MAT) ! out of date material recovery, see FLT_MAT !! * * * * * * * * * * * * * * * * * * * * * * * * * * !! EXTRACT REAL MATERIAL PROPERTIES OF MATERIAL !! NUM FROM MISCELLANEOUS REAL SYSTEM PROPERTIES !! * * * * * * * * * * * * * * * * * * * * * * * * * * !Use System_Constants !IMPLICIT NONE ! INTEGER, INTENT(IN) :: NUM, N_PR_MAT ! REAL(DP), INTENT(IN) :: FLT_MISC (MISC_FL) ! REAL(DP), INTENT(OUT) :: PRT_MAT (N_PR_MAT) ! INTEGER :: I, ISTART, N_MAX ! !! NUM = MATERIAL NUMBER !! N_MAX = MAXIMUM ALLOWABLE MATERIAL NUMBER !! N_PR_MAT = NUMBER OF REAL PROPERTIES PER MATERIAL !! MISC_FL = NUMBER OF MISC REAL SYSTEM PROPERTIES !! FLT_MISC = SYSTEM STORAGE FOR MISC REAL PROP !! PRT_MAT = REAL PROPERTY ARRAY FOR MATERIAL NUM, ! !! PROPERTIES STORED IN FLT_MISC IN ORDER OF MATERIAL NUMBER ! IF ( N_PR_MAT < 1 .OR. N_PR_MAT > MISC_FL ) & ! STOP 'BAD N_PR_MAT, IN GET_MATERIAL_PROPS' ! N_MAX = MISC_FL / N_PR_MAT ! Max material number ! IF ( NUM > N_MAX ) STOP 'DIMENSIONS EXCEEDED, GET_MATERIAL_PROPS' ! IF ( N_MAX * N_PR_MAT /= MISC_FL ) PRINT *, & ! 'WARNING: MISC_FL > MATERIAL PROPERTY STORAGE' ! N_WARN = N_WARN + 1 ! INCREMENT WARNING ! ISTART = N_PR_MAT * (NUM - 1) ! DO I = 1, N_PR_MAT ! PRT_MAT (I) = FLT_MISC (ISTART + I) ! END DO ! END SUBROUTINE GET_MATERIAL_PROPS ! SUBROUTINE LIST_ALL_MATERIALS (N_PR_MAT, FLT_MISC) ! out of date method ! !b SUBROUTINE LIST_ALL_MATERIALS (N_PR_MAT, FLT_MISC, PRT_MAT) !! * * * * * * * * * * * * * * * * * * * * * * * * * !! LIST REAL PROPERTIES BY MATERIAL NUMBER !! * * * * * * * * * * * * * * * * * * * * * * * * * !Use System_Constants !IMPLICIT NONE ! INTEGER, INTENT(IN) :: N_PR_MAT ! REAL(DP), INTENT(IN) :: FLT_MISC (MISC_FL) ! !b REAL(DP), INTENT(IN) :: PRT_MAT (MISC_FL) ! REAL(DP) :: PRT_MAT (N_PR_MAT) ! INTEGER :: I, J, N_MAX ! !! N_PR_MAT = NUMBER OF REAL PROPERTIES PER MATERIAL !! <= MISC_FL !! MISC_FL = NUMBER OF MISC REAL SYSTEM PROP !! FLT_MISC = SYSTEM STORAGE OF MISC REAL PROP !! PRT_MAT = REAL PROP ARRAY FOR MATERIAL NUM !! N_MAX = MAXIMUM ALLOWABLE MATERIAL NUMBER ! ! IF ( N_PR_MAT < 1 .OR. N_PR_MAT > MISC_FL) THEN ! WRITE (N_PRT, * ) 'N_PR_MAT, MISC_FL ', N_PR_MAT, MISC_FL ! STOP 'BAD N_PR_MAT, IN LIST_ALL_MATERIALS' ! ENDIF ! N_MAX = MISC_FL / N_PR_MAT ! WRITE (N_PRT, 50) N_PR_MAT ! 50 FORMAT (/, '** LIST OF REAL MATERIAL PROPERTIES **',/, & ! 'MATERIAL,',I4,' REAL PROPERTIES') ! DO I = 1, N_MAX ! CALL GET_MATERIAL_PROPS (I, N_PR_MAT, FLT_MISC, PRT_MAT) ! WRITE (N_PRT, 5010) I, (PRT_MAT (J), J = 1, N_PR_MAT) ! 5010 FORMAT ( I5, (10(1PE12.4)) ) ! END DO ! END SUBROUTINE LIST_ALL_MATERIALS SUBROUTINE PT_CODE (JPT, KODE, KODES) ! * * * * * * * * * * * * * * * * * * * * * * * * * * ! EXTRACT B.C. INDICATORS AT NODE NUMBER JPT ! * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants IMPLICIT NONE INTEGER, INTENT(IN) :: JPT, KODE INTEGER, INTENT(OUT) :: KODES (N_G_DOF) INTEGER :: I, II, IK, INEW, IOLD, ISUM ! JPT = NODE NO. ! N_G_DOF = NO. PARAMETERS PER NODE ! KODE = (N_G_DOF) DIGIT INTEGER CONTAINING BC INDICATORS ! KODES = VECTOR CONTAINING N_G_DOF INTEGER CODES (0 OR I) ! 0 IMPLIES NO B. C. ! J IMPLIES A B. C. OF TYPE J IOLD = KODE ISUM = 0 DO I = 1, N_G_DOF II = N_G_DOF + 1 - I INEW = IOLD / 10 IK = IOLD - INEW * 10 ISUM = ISUM + IK * 10** (I - 1) IOLD = INEW KODES (II) = IK IF ( IK > MAX_TYP ) THEN PRINT *, 'INVALID BC CODE, NODE = ', JPT, ', CODE = ', KODE PRINT *, 'INVALID BC TYPE, NODE = ', JPT, ', PARM = ', I PRINT *, 'Set keyword MAX_TYP to ', IK STOP 'INVALID BC TYPE, PT_CODE' END IF END DO ! WAS DATA RIGHT JUSTIFIED? IF ( KODE > ISUM ) THEN WRITE (N_BUG, * ) 'WARNING: BC NOT RIGHT JUSTIFIED AT NODE ', JPT N_WARN = N_WARN + 1 ! INCREMENT WARNING END IF END SUBROUTINE PT_CODE SUBROUTINE COUNT_L_AT_NODE (N_ELEMS, NOD_PER_EL, MAX_NP, & NODES, L_TO_N_SUM) ! * * * * * * * * * * * * * * * * * * * * * * * * * * ! COUNT NUMBER OF ELEMENTS ADJACENT TO EACH NODE ! (TO SIZE ELEM ADJACENT TO ELEM LIST) ! * * * * * * * * * * * * * * * * * * * * * * * * * * IMPLICIT NONE INTEGER, INTENT(IN) :: N_ELEMS, NOD_PER_EL, MAX_NP INTEGER, INTENT(IN) :: NODES (N_ELEMS, NOD_PER_EL) INTEGER, INTENT(OUT) :: L_TO_N_SUM (MAX_NP) INTEGER :: ELEM_NODES (NOD_PER_EL) INTEGER :: IE, IN, N_TEMP ! ELEM_NODES = INCIDENCES ARRAY FOR A SINGLE ELEMENT ! L_TO_N_SUM (I) = NUMBER OF ELEM NEIGHBORS OF NODE I ! MAX_NP = TOTAL NUMBER OF NODES ! N_ELEMS = TOTAL NUMBER OF ELEMENTS ! NOD_PER_EL = NUMBER OF NODES PER ELEMENT ! NODES = SYSTEM ARRAY OF INCIDENCES OF ALL ELEMENTS L_TO_N_SUM = 0 ! INITIALIZE DO IE = 1, N_ELEMS ! LOOP OVER ALL ELEMENTS ! EXTRACT INCIDENCES LIST FOR ELEMENT IE CALL ELEMENT_NODES (IE, NOD_PER_EL, NODES, ELEM_NODES) DO IN = 1, NOD_PER_EL ! loop over each node N_TEMP = ELEM_NODES (IN) IF ( N_TEMP < 1 ) CYCLE ! to a real node L_TO_N_SUM (N_TEMP) = L_TO_N_SUM (N_TEMP) + 1 END DO ! over IN END DO ! over elements END SUBROUTINE COUNT_L_AT_NODE SUBROUTINE FORM_L_ADJACENT_NODES (N_ELEMS, NOD_PER_EL, MAX_NP, & NODES, NEIGH_N, L_TO_N_NEIGH) ! * * * * * * * * * * * * * * * * * * * * * * * * * * ! TABULATE ELEMENTS ADJACENT TO EACH NODE ! (TO SIZE ELEM ADJACENT TO ELEM LIST) ! * * * * * * * * * * * * * * * * * * * * * * * * * * IMPLICIT NONE INTEGER, INTENT(IN) :: N_ELEMS, NOD_PER_EL, MAX_NP, NEIGH_N INTEGER, INTENT(IN) :: NODES (N_ELEMS, NOD_PER_EL) INTEGER, INTENT(OUT) :: L_TO_N_NEIGH (NEIGH_N, MAX_NP) INTEGER :: ELEM_NODES (NOD_PER_EL) ! scratch area INTEGER :: COUNT (MAX_NP) ! scratch area INTEGER :: IE, IN, N_TEMP ! ELEM_NODES = INCIDENCES ARRAY FOR A SINGLE ELEMENT ! L_TO_N_SUM (I) = NUMBER OF ELEM NEIGHBORS OF NODE I ! MAX_NP = TOTAL NUMBER OF NODES ! N_ELEMS = TOTAL NUMBER OF ELEMENTS ! NOD_PER_EL = NUMBER OF NODES PER ELEMENT ! NODES = SYSTEM ARRAY OF INCIDENCES OF ALL ELEMENTS L_TO_N_NEIGH = 0 ! INITIALIZE COUNT = 0 ! INITIALIZE DO IE = 1, N_ELEMS ! LOOP OVER ALL ELEMENTS ! EXTRACT INCIDENCES LIST FOR ELEMENT IE CALL ELEMENT_NODES (IE, NOD_PER_EL, NODES, ELEM_NODES) DO IN = 1, NOD_PER_EL ! loop over each node N_TEMP = ELEM_NODES (IN) IF ( N_TEMP < 1 ) CYCLE ! to a real node COUNT (N_TEMP) = COUNT (N_TEMP) + 1 IF ( COUNT (N_TEMP) > NEIGH_N ) & STOP 'INCONSISTENT NEIGH_N; FORM_L_ADJACENT_NODES' L_TO_N_NEIGH (COUNT (N_TEMP), N_TEMP) = IE END DO ! over IN END DO ! over elements END SUBROUTINE FORM_L_ADJACENT_NODES !SUBROUTINE COUNT_L_AT_ELEM (N_ELEMS, NOD_PER_EL, MAX_NP, & ! NODES, L_TO_N_SUM, L_TO_L_SUM) !! * * * * * * * * * * * * * * * * * * * * * * * * * * !! -------This is incorrect except for quadratic elements------- !! COUNT NUMBER OF ELEMENTS ADJACENT TO OTHER ELEMENTS !! (SEE ALTERNATE: COUNT_ELEMS_AT_ELEM) !! * * * * * * * * * * * * * * * * * * * * * * * * * * ! IMPLICIT NONE ! INTEGER, INTENT(IN) :: N_ELEMS, NOD_PER_EL, MAX_NP ! INTEGER, INTENT(IN) :: NODES (N_ELEMS, NOD_PER_EL) ! INTEGER, INTENT(OUT) :: L_TO_N_SUM (MAX_NP) ! INTEGER, INTENT(OUT) :: L_TO_L_SUM (N_ELEMS) ! ! INTEGER :: ELEM_NODES (NOD_PER_EL) ! INTEGER :: FOUND, IE, IN, N_TEST ! !! ELEM_NODES = INCIDENCES ARRAY FOR A SINGLE ELEMENT !! L_TO_L_SUM (I) = NUMBER OF ELEM NEIGHBORS OF ELEMENT I !! L_TO_N_SUM (I) = NUMBER OF ELEM NEIGHBORS OF NODE I !! MAX_NP = TOTAL NUMBER OF NODES !! N_ELEMS = TOTAL NUMBER OF ELEMENTS !! NOD_PER_EL = NUMBER OF NODES PER ELEMENT !! NODES = SYSTEM ARRAY OF INCIDENCES OF ALL ELEMENTS ! ! DO IE = 1, N_ELEMS ! LOOP OVER ALL ELEMENTS !if (ie ==1 ) stop'COUNT_L_AT_ELEM is incorrect' ! FOUND = 0 ! !! EXTRACT INCIDENCES LIST FOR ELEMENT IE ! CALL ELEMENT_NODES (IE, NOD_PER_EL, NODES, ELEM_NODES) ! ! DO IN = 1, NOD_PER_EL ! LOOP OVER NODES OF ELEMENT ! N_TEST = ELEM_NODES (IN) ! IF ( N_TEST > 0 ) FOUND = FOUND + L_TO_N_SUM (N_TEST) ! END DO ! over in ! !! REMOVE ELEMENTS COUNTED TWICE ! FOUND = FOUND - 2 * COUNT ( ELEM_NODES > 0 ) ! L_TO_L_SUM (IE) = FOUND ! END DO ! over all elements !END SUBROUTINE COUNT_L_AT_ELEM SUBROUTINE COUNT_ELEMS_AT_ELEM (N_ELEMS, NOD_PER_EL, MAX_NP, & L_FIRST, L_LAST, NODES, NEEDS, & L_TO_L_SUM, N_WARN) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! COUNT NUMBER OF ELEMENTS ADJACENT TO OTHER ELEMENTS ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * IMPLICIT NONE INTEGER, INTENT(IN) :: N_ELEMS, NOD_PER_EL, MAX_NP, NEEDS INTEGER, INTENT(IN) :: L_FIRST (MAX_NP), L_LAST (MAX_NP) INTEGER, INTENT(IN) :: NODES (N_ELEMS, NOD_PER_EL) INTEGER, INTENT(OUT) :: L_TO_L_SUM (N_ELEMS) INTEGER, INTENT(INOUT) :: N_WARN INTEGER :: ELEM_NODES (NOD_PER_EL) INTEGER :: NEIG_NODES (NOD_PER_EL) INTEGER :: FOUND, IE, IN, L_TEST, L_START, L_STOP, N_TEST, NULLS INTEGER :: NEED, KOUNT ! ELEM_NODES = INCIDENCES ARRAY FOR A SINGLE ELEMENT ! KOUNT = NUMBER OF COMMON NODES ! L_FIRST (I) = ELEMENT WHERE NODE I FIRST APPEARS ! L_LAST (I) = ELEMENT WHERE NODE I LAST APPEARS ! L_TO_L_SUM (I) = NUMBER OF ELEM NEIGHBORS OF ELEMENT I ! MAX_NP = TOTAL NUMBER OF NODES ! NEEDS = NUMBER OF COMMON NODES TO BE A NEIGHBOR ! N_ELEMS = TOTAL NUMBER OF ELEMENTS ! NOD_PER_EL = NUMBER OF NODES PER ELEMENT ! NODES = SYSTEM ARRAY OF INCIDENCES OF ALL ELEMENTS NEED = MAX (1, NEEDS) ! INITIALIZE L_TO_L_SUM = 0 MAIN : DO IE = 1, N_ELEMS ! LOOP OVER ALL ELEMENTS FOUND = 0 ! EXTRACT INCIDENCES LIST FOR ELEMENT IE CALL ELEMENT_NODES (IE, NOD_PER_EL, NODES, ELEM_NODES) ! ESTABLISH RANGE OF POSSIBLE ELEMENT NEIGHBORS L_START = N_ELEMS ; L_STOP = 0 DO IN = 1, NOD_PER_EL N_TEST = ELEM_NODES (IN) !b IF ( N_TEST < 1 ) CYCLE! to a real node !b L_START = MIN (L_START, L_FIRST (N_TEST) ) !b L_STOP = MAX (L_STOP, L_LAST (N_TEST) ) !b END DO ! LOOP OVER POSSIBLE ELEMENT NEIGHBORS IF ( L_START <= L_STOP) THEN RANGE : DO L_TEST = L_START, L_STOP IF ( L_TEST /= IE) THEN KOUNT = 0 ! NO COMMON NODES ! LOOP OVER INCIDENCES OF POSSIBLE ELEMENT NEIGHBOR CALL ELEMENT_NODES (L_TEST, NOD_PER_EL, NODES, NEIG_NODES) LOCAL : DO IN = 1, NOD_PER_EL N_TEST = NEIG_NODES (IN) IF ( N_TEST < 1 .OR. N_TEST > MAX_NP ) THEN PRINT *, 'COUNT_ELEMS_AT_ELEM: WARNING, ', & 'SKIPPING INVALID NODE ', N_TEST, & ' OF ELEMENT ', L_TEST PRINT *, NEIG_NODES N_WARN = N_WARN + 1 ! INCREMENT WARNING CYCLE LOCAL ! to a real node END IF ! IMPOSSIBLE NODE IF ( L_FIRST (N_TEST) > IE ) CYCLE LOCAL ! to next node IF ( L_LAST (N_TEST) < IE ) CYCLE LOCAL ! to next node ! COMPARE WITH INCIDENCES OF ELEMENT IE IF ( ANY ( ELEM_NODES == N_TEST ) ) THEN KOUNT = KOUNT + 1 IF ( KOUNT == NEED ) THEN ! IS A NEIGHBOR FOUND = FOUND + 1 EXIT LOCAL ! this L_TEST element search loop END IF ! NUMBER NEEDED END IF END DO LOCAL ! over in END IF END DO RANGE ! over candidate element L_TEST END IF ! a possible candidate L_TO_L_SUM (IE) = FOUND END DO MAIN ! over all elements PRINT *, 'NOTE: MAXIMUM NUMBER OF ELEMENT NEIGHBORS = ', & MAXVAL ( L_TO_L_SUM ) NULLS = COUNT ( L_TO_L_SUM == 0 ) ! CHECK DATA IF ( NULLS > 0 ) THEN PRINT *, 'WARNING, ', NULLS, ' ELEMENTS HAVE NO ELEMENT NEIGHBORS' N_WARN = N_WARN + 1 ! INCREMENT WARNING END IF END SUBROUTINE COUNT_ELEMS_AT_ELEM SUBROUTINE FIRST_LAST_L_AT_PT (NODES, L_FIRST, L_LAST, I_WARN) !* * * * * * * * * * * * * * * * * * * * * * * * * * * * ! FIND ELEMENT OF FIRST & LAST APPEARANCE OF EACH NODE !* * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants Use Elem_type_data Use Interface_header IMPLICIT NONE INTEGER, INTENT(IN) :: I_WARN INTEGER, INTENT(IN) :: NODES (L_S_TOT, NOD_PER_EL) INTEGER, INTENT(OUT) :: L_FIRST (MAX_NP), L_LAST (MAX_NP) INTEGER, ALLOCATABLE :: EL_NODES (:) INTEGER :: I, IE, L, LT ! EL_NODES = INCIDENCES ARRAY OF A SINGLE ELEMENT ! I_WARN > 0 WARN OF NODES WITH NO ELEMENTS ! L_FIRST(I) = ELEMENT OF FIRST APPEARANCE OF NODE I ! L_LAST(I) = ELEMENT OF LAST APPEARANCE OF NODE I ! LT_N = MAXIMUM NUMBER OF NODES FOR ELEMENT TYPE ! L_S_TOT = TOTAL GENERALIZED ELEMENTS, N_ELEMS + N_MIXED + N_F_SEG ! MAX_NP = TOTAL NUMBER OF NODES IN SYSTEM ! NODES = SYSTEN ARRAY OF ALL ELEMENT INCIDENCES ! NOD_PER_EL = NUMBER OF NODES PER ELEMENT ! N_ELEMS = NUMBER OF ELEMENTS IN THE SYSTEM ! N_L_TYPE = NUMBER OF DIFFERENT ELEMENT TYPES USED L_FIRST = 0 ; L_LAST = 0 ! initialize LT = 1 ! INITIALIZE ELEMENT TYPE DO IE = 1, N_ELEMS !--> 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 GET_ELEM_TYPE_DATA (LT) !b SET_ELEM_TYPE_INFO (LT) LAST_LT = LT !--> Allocate Arrays for this type IF ( IE > 1 ) DEALLOCATE ( EL_NODES ) !bu ALLOCATE ( EL_NODES (LT_N) ) !bu END IF ! a new element type ! EXTRACT ELEMENT'S NODES EL_NODES = GET_ELEM_NODES (IE, LT_N, NODES) ! SCAN THE NODES DO I = 1, LT_N L = EL_NODES (I) IF ( L < 1 ) CYCLE ! to next node IF ( L_FIRST (L) == 0 ) L_FIRST (L) = IE L_LAST (L) = IE END DO END DO DEALLOCATE ( EL_NODES ) ! WARN OF NODES WITH NO ELEMENT CONNECTIONS IF ( I_WARN == 0 .OR. COUNT( L_FIRST == 0 ) == 0 ) RETURN DO I = 1, MAX_NP IF ( L_FIRST (I) > 0 ) CYCLE ! to next node PRINT *, ' ' WRITE (6, 5) I ; 5 FORMAT ('WARNING, NODE ', I6, & & ' DOES NOT OCCUR IN THE ELEMENT CONNECTION LIST') N_WARN = N_WARN + 1 ! INCREMENT WARNING END DO END SUBROUTINE FIRST_LAST_L_AT_PT SUBROUTINE LIST_FIRST_LAST (L_FIRST, L_LAST) !* * * * * * * * * * * * * * * * * * * * * * * * * * * * ! LIST ELEMENT OF FIRST & LAST APPEARANCE OF EACH NODE !* * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants IMPLICIT NONE INTEGER, INTENT(IN) :: L_FIRST (MAX_NP), L_LAST (MAX_NP) INTEGER :: I ! L_FIRST(I) = ELEMENT OF FIRST APPEARANCE OF NODE I ! L_LAST(I) = ELEMENT OF LAST APPEARANCE OF NODE I ! MAX_NP = TOTAL NUMBER OF NODES IN SYSTEM WRITE (N_PRT, *) ' ' WRITE (N_PRT, *) '*** ELEMENTS AT NODE LIST ***' WRITE (N_PRT, *) 'NODE FIRST_ELEM LAST_ELEM' DO I = 1, MAX_NP WRITE (N_PRT, '(I6, 5X, I6, 5X, I6)') & I, L_FIRST (I), L_LAST (I) END DO WRITE (N_PRT, *) ' ' END SUBROUTINE LIST_FIRST_LAST SUBROUTINE LIST_L_AT_NODE_INFO (L_TO_N_SUM, L_FIRST, L_LAST) !* * * * * * * * * * * * * * * * * * * * * * * * * * * * ! LIST ELEMENT OF FIRST & LAST APPEARANCE OF EACH NODE ! AND THE NUMBER OF ELEMENTS ATTACHED !* * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants IMPLICIT NONE INTEGER, INTENT(IN) :: L_FIRST (MAX_NP), L_LAST (MAX_NP) INTEGER, INTENT(IN) :: L_TO_N_SUM (MAX_NP) INTEGER :: I ! L_FIRST(I) = ELEMENT OF FIRST APPEARANCE OF NODE I ! L_LAST(I) = ELEMENT OF LAST APPEARANCE OF NODE I ! L_TO_N_SUM (I) = NUMBER OF ELEM NEIGHBORS OF NODE I ! MAX_NP = TOTAL NUMBER OF NODES IN SYSTEM WRITE (N_PRT, *) ' ' WRITE (N_PRT, *) '*** ELEMENTS AT NODE SUMMARY ***' WRITE (N_PRT, *) ' NODE TOTAL_ELEM FIRST_ELEM LAST_ELEM' DO I = 1, MAX_NP WRITE (N_PRT, '(I6, 1X, I6, 5X, I6, 5X, I6)') & I, L_TO_N_SUM (I), L_FIRST (I), L_LAST (I) END DO WRITE (N_PRT, *) ' ' END SUBROUTINE LIST_L_AT_NODE_INFO SUBROUTINE LIST_ELEMENTS_AT_EL (L_TO_L_NEIGH) !* * * * * * * * * * * * * * * * * * * * * * * * * * * * ! LIST ELEMENT NEIGHBORS OF EACH ELEMENT !* * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants ! for NEIGH_L, N_ELEMS IMPLICIT NONE INTEGER, INTENT(IN) :: L_TO_L_NEIGH (NEIGH_L, N_ELEMS) INTEGER :: IE ! L_TO_L_NEIGH (J, I) = ELEMENT NEIGHBOR J OF ELEMENT I ! NEIGH_L = MAX NUMBER OF ELEM NEIGHBORS AT ANY ELEM ! N_ELEMS = MAXIMUM NUMBER OF ELEMENTS IN SYSTEM WRITE (N_PRT, *) ' ' WRITE (N_PRT, *) '*** NEIGHBORING ELEMENTS AT EACH ELEMENT ***' WRITE (N_PRT, *) 'ELEMENT, ADJACENT ELEMENT LIST' DO IE = 1, N_ELEMS WRITE (N_PRT, 5) IE, L_TO_L_NEIGH (1:NEIGH_L, IE) 5 FORMAT (I6, ',', 10I6, (/, 7X, 10I6) ) END DO WRITE (N_PRT, *) ' ' END SUBROUTINE LIST_ELEMENTS_AT_EL SUBROUTINE LIST_ELEMENTS_AT_PT (L_TO_N_NEIGH) !* * * * * * * * * * * * * * * * * * * * * * * * * * * * ! LIST ELEMENT NEIGHBORS OF EACH NODE !* * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants ! for NEIGH_N, MAX_NP IMPLICIT NONE INTEGER, INTENT(IN) :: L_TO_N_NEIGH (NEIGH_N, MAX_NP) INTEGER :: IP ! L_TO_N_NEIGH (J, I) = ELEMENT NEIGHBOR J OF NODE PT I ! NEIGH_N = MAXIMUM NUMBER OF ELEM NEIGHBORS AT ANY ELEM ! MAX_NP = MAXIMUM NUMBER OF NODE PTS IN SYSTEM WRITE (N_PRT, *) ' ' WRITE (N_PRT, *) '*** NEIGHBORING ELEMENTS AT EACH NODE ***' WRITE (N_PRT, *) ' NODE, ADJACENT ELEMENT LIST' DO IP = 1, MAX_NP WRITE (N_PRT, 5) IP, L_TO_N_NEIGH (1:NEIGH_N, IP) 5 FORMAT (I8, ', ', 10I6, (/, 12X, 10I6) ) END DO WRITE (N_PRT, *) ' ' END SUBROUTINE LIST_ELEMENTS_AT_PT SUBROUTINE FORM_ELEMS_AT_EL (N_ELEMS, NOD_PER_EL, MAX_NP, & L_FIRST, L_LAST, NODES, N_SPACE, & L_TO_L_SUM, L_TO_L_NEIGH, & NEIGH_L, NEEDS, ON_BOUNDARY, U_ON_B) ! * * * * * * * * * * * * * * * * * * * * * * * * * * ! FORM LIST OF ELEMENTS ADJACENT TO OTHER ELEMENTS ! * * * * * * * * * * * * * * * * * * * * * * * * * * IMPLICIT NONE INTEGER, INTENT(IN) :: N_ELEMS, NOD_PER_EL, MAX_NP, NEIGH_L INTEGER, INTENT(IN) :: L_FIRST (MAX_NP), L_LAST (MAX_NP) INTEGER, INTENT(IN) :: NODES (N_ELEMS, NOD_PER_EL) INTEGER, INTENT(IN) :: L_TO_L_SUM (N_ELEMS) INTEGER, INTENT(IN) :: N_SPACE, NEEDS, U_ON_B INTEGER, INTENT(OUT) :: L_TO_L_NEIGH (NEIGH_L, N_ELEMS) LOGICAL, INTENT(INOUT) :: ON_BOUNDARY (N_ELEMS) INTEGER :: ELEM_NODES (NOD_PER_EL) INTEGER :: NEIG_NODES (NOD_PER_EL) INTEGER :: IE, IN, L_TEST, L_START, L_STOP, N_TEST INTEGER :: FOUND, NEXT ! add NEED for face or edge sub-sets INTEGER :: SUM_L_TO_L INTEGER :: IO_1, KOUNT, NEED, N_FACES, WHERE !b ! ON_BOUNDARY = TRUE IF ELEMENT HAS A FACE ON THE BOUNDARY ! ELEM_NODES = INCIDENCES ARRAY FOR A SINGLE ELEMENT ! FOUND = CURRENT NUMBER OF LOCAL NEIGHBORS ! KOUNT = CURRENT NUMBER OF COMMON NODES ! L_FIRST (I) = ELEMENT WHERE NODE I FIRST APPEARS ! L_LAST (I) = ELEMENT WHERE NODE I LAST APPEARS ! L_TO_L_NEIGH (J, I) = ELEM NEIGHBOR J OF ELEMENT I ! L_TO_L_SUM (I) = NUMBER OF ELEM NEIGHBORS OF ELEMENT I ! NEEDS = NUMBER OF COMMON NODES TO BE A NEIGHBOR ! NEIGH_L = MAXIMUM NUMBER OF NEIGHBORS AT AN ELEMENT ! MAX_NP = TOTAL NUMBER OF NODES ! N_ELEMS = TOTAL NUMBER OF ELEMENTS ! NOD_PER_EL = NUMBER OF NODES PER ELEMENT ! NODES = SYSTEM ARRAY OF INCIDENCES OF ALL ELEMENTS ! WHERE = LOCATION TO INSERT NEIGHBOR, <= MAX_FACES NEED = MAX (1, NEEDS) ! INITIALIZE L_TO_L_NEIGH = 0 ! INITIALIZE MAIN : DO IE = 1, N_ELEMS ! LOOP OVER ALL ELEMENTS SUM_L_TO_L = L_TO_L_SUM (IE) ! MAX POSSIBLE NEIGHBORS FOUND = COUNT ( L_TO_L_NEIGH (:, IE) > 0 ) ! PREVIOUSLY FOUND IF ( FOUND == SUM_L_TO_L ) CYCLE MAIN ! ALL HAVE BEEN FOUND ! EXTRACT INCIDENCES LIST FOR ELEMENT IE CALL ELEMENT_NODES (IE, NOD_PER_EL, NODES, ELEM_NODES) ! ESTABLISH RANGE OF POSSIBLE ELEMENT NEIGHBORS L_START = N_ELEMS + 1 ; L_STOP = 0 DO IN = 1, NOD_PER_EL L_TEST = ELEM_NODES (IN) !b IF ( L_TEST < 1 ) CYCLE ! TO NEXT NODE !b L_START = MIN (L_START, L_FIRST (L_TEST)) !b L_STOP = MAX (L_STOP, L_LAST (L_TEST)) !b !b L_START = MIN (L_START, L_FIRST (ELEM_NODES (IN)) ) !b L_STOP = MAX (L_STOP, L_LAST (ELEM_NODES (IN)) ) END DO L_START = MAX (L_START, IE+1) ! SEARCH ABOVE IE ONLY ! LOOP OVER POSSIBLE ELEMENT NEIGHBORS IF ( L_START <= L_STOP) THEN RANGE : DO L_TEST = L_START, L_STOP KOUNT = 0 ! NO COMMON NODES ! EXTRACT NODES OF L_TEST CALL ELEMENT_NODES (L_TEST, NOD_PER_EL, NODES, NEIG_NODES) ! LOOP OVER INCIDENCES OF POSSIBLE ELEMENT NEIGHBOR LOCAL : DO IN = 1, NOD_PER_EL N_TEST = NEIG_NODES (IN) IF ( N_TEST < 1 ) CYCLE ! to a real node IF ( L_FIRST (N_TEST) > IE ) CYCLE LOCAL ! to next node IF ( L_LAST (N_TEST) < IE ) CYCLE LOCAL ! to next node ! COMPARE WITH INCIDENCES OF ELEMENT IE IF ( ANY ( ELEM_NODES == N_TEST ) ) THEN KOUNT = KOUNT + 1 ! SHARED NODE COUNT IF ( KOUNT == NEED ) THEN ! NEIGHBOR PAIR FOUND FOUND = FOUND + 1 ! INSERT THE PAIR ! NOTE: THIS INSERT IS NOT ORDERED. IF A SPECIFIC FACE ! ORDER IS REQUIRED DETERMINE WHERE HERE WHERE = FOUND ! OR ORDER THE CURRENT FACE L_TO_L_NEIGH (WHERE, IE) = L_TEST ! ONE OF TWO NEXT = COUNT ( L_TO_L_NEIGH(:, L_TEST) > 0 ) !b /= ?? WHERE = NEXT+1 ! OR ORDER THE NEIGHBOR FACE IF ( L_TO_L_SUM (L_TEST) > NEXT ) & L_TO_L_NEIGH (NEXT+1, L_TEST) = IE ! TWO OF TWO IF ( SUM_L_TO_L == FOUND ) CYCLE MAIN ! ALL FOUND CYCLE RANGE ! this L_TEST element search loop END IF ! NUMBER NEEDED END IF ! SHARE AT LEAST ONE COMMON NODE END DO LOCAL ! over N_TEST END DO RANGE ! over candidate element L_TEST END IF ! a possible candidate END DO MAIN! over all elements IF ( NEED >= N_SPACE ) THEN ! DATA ARE EDGE OR FACE NEIGHBORS !b call at(2698) ON_BOUNDARY = .FALSE. ! SAVE THE ELEMENT NUMBERS THAT FACE THE BOUNDARY DO IE = 1, N_ELEMS CALL GET_LT_FACES (IE, N_FACES) IF ( N_FACES > 0 ) THEN ! MIGHT BE ON THE BOUNDARY IF ( NEIGH_L < N_FACES ) THEN PRINT *,'NOTE: REDUCING FACES FROM ', N_FACES, & & ' TO ', NEIGH_L N_FACES = NEIGH_L ! Possibly valid mesh END IF ! validation IF ( ANY (L_TO_L_NEIGH (1:N_FACES, IE) == 0) ) THEN ! YES ON_BOUNDARY (IE) = .TRUE. END IF ! ON BOUNDARY END IF ! POSSIBLE ELEMENT END DO ! OVER ELEMENTS END IF ! SEARCH OF FACING NEIGHBORS IF ( ANY ( ON_BOUNDARY ) ) THEN ! SAVE ELEMENTS ON BOUNDARY OPEN (U_ON_B, FILE='boundary_face_elems.tmp', & ACTION='WRITE', STATUS='REPLACE', IOSTAT = IO_1) IF ( IO_1 > 0 ) PRINT *, 'OPEN FAILED FOR boundary_face_elems.tmp' DO IE = 1, N_ELEMS IF ( ON_BOUNDARY (IE) ) WRITE (U_ON_B, *) IE END DO PRINT *,'NOTE: Saved ', COUNT ( ON_BOUNDARY ), & ' elements on the boundary.' END IF ! SAVE ELEMENTS ON BOUNDARY END SUBROUTINE FORM_ELEMS_AT_EL !!SUBROUTINE FORM_ELEMS_AT_EL_OLD (N_ELEMS, NOD_PER_EL, MAX_NP, & ! L_FIRST, L_LAST, NODES, & ! L_TO_L_SUM, L_TO_L_NEIGH, NEIGH_L) !! * * * * * * * * * * * * * * * * * * * * * * * * * * !! FORM LIST OF ELEMENTS ADJACENT TO OTHER ELEMENTS !! * * * * * * * * * * * * * * * * * * * * * * * * * * ! IMPLICIT NONE ! INTEGER, INTENT(IN) :: N_ELEMS, NOD_PER_EL, MAX_NP, NEIGH_L ! INTEGER, INTENT(IN) :: L_FIRST (MAX_NP), L_LAST (MAX_NP) ! INTEGER, INTENT(IN) :: NODES (N_ELEMS, NOD_PER_EL) ! INTEGER, INTENT(IN) :: L_TO_L_SUM (N_ELEMS) ! INTEGER, INTENT(OUT) :: L_TO_L_NEIGH (NEIGH_L, N_ELEMS) ! ! INTEGER :: ELEM_NODES (NOD_PER_EL) ! INTEGER :: NEIG_NODES (NOD_PER_EL) ! INTEGER :: IE, IN, L_TEST, L_START, L_STOP, N_TEST ! INTEGER :: FOUND ! !! ELEM_NODES = INCIDENCES ARRAY FOR A SINGLE ELEMENT !! FOUND = CURRENT NUMBER OF LOCAL NEIGHBORS !! L_FIRST (I) = ELEMENT WHERE NODE I FIRST APPEARS !! L_LAST (I) = ELEMENT WHERE NODE I LAST APPEARS !! L_TO_L_NEIGH (J, I) = ELEM NEIGHBOR J OF ELEMENT I !! L_TO_L_SUM (I) = NUMBER OF ELEM NEIGHBORS OF ELEMENT I !! NEIGH_L = MAXIMUM NUMBER OF NEIGHBORS AT AN ELEMENT !! MAX_NP = TOTAL NUMBER OF NODES !! N_ELEMS = TOTAL NUMBER OF ELEMENTS !! NOD_PER_EL = NUMBER OF NODES PER ELEMENT !! NODES = SYSTEM ARRAY OF INCIDENCES OF ALL ELEMENTS ! ! L_TO_L_NEIGH = 0 ! INITIALIZE ! DO IE = 1, N_ELEMS ! LOOP OVER ALL ELEMENTS ! FOUND = 0 ! !! EXTRACT INCIDENCES LIST FOR ELEMENT IE ! CALL ELEMENT_NODES (IE, NOD_PER_EL, NODES, ELEM_NODES) ! !! ESTABLISH RANGE OF POSSIBLE ELEMENT NEIGHBORS ! L_START = HUGE (L_START) ; L_STOP = 0 ! DO IN = 1, NOD_PER_EL !! XXX will get zero for mixed element types XXXX ! IF ( L_START > L_FIRST (ELEM_NODES (IN) ) ) & ! L_START = L_FIRST (ELEM_NODES (IN) ) ! IF ( L_STOP < L_LAST (ELEM_NODES (IN) ) ) & ! L_STOP = L_LAST (ELEM_NODES (IN) ) ! END DO ! !! LOOP OVER POSSIBLE ELEMENT NEIGHBORS ! IF ( L_START <= L_STOP) THEN ! DO L_TEST = L_START, L_STOP ! IF ( L_TEST /= IE) THEN ! !! LOOP OVER INCIDENCES OF POSSIBLE ELEMENT NEIGHBOR ! CALL ELEMENT_NODES (L_TEST, NOD_PER_EL, NODES, NEIG_NODES) ! DO IN = 1, NOD_PER_EL ! N_TEST = NEIG_NODES (IN) ! IF ( N_TEST < 1 ) CYCLE ! to a real node ! IF ( L_FIRST (N_TEST) > IE ) CYCLE ! to next node ! IF ( L_LAST (N_TEST) < IE ) CYCLE ! to next node ! !! COMPARE WITH INCIDENCES OF ELEMENT IE ! IF ( ANY ( ELEM_NODES == N_TEST ) ) THEN ! FOUND = FOUND + 1 ! L_TO_L_NEIGH (FOUND, IE) = L_TEST ! EXIT ! this L_TEST element search loop ! END IF ! END DO ! over N_TEST ! END IF ! not itself ! IF ( L_TO_L_SUM (IE) == FOUND ) EXIT ! have last candidate ! END DO ! over candidate element L_TEST ! END IF ! a possible candidate ! END DO ! over all elements !END SUBROUTINE FORM_ELEMS_AT_EL_OLD SUBROUTINE SAVE_QP_LOOP_COUNT (ELEM_QP) ! .............................................................. ! SAVE NUMBER OF QUADRATURE POINTS THAT AN ELEMENT IS ! USING (FOR SCP RECOVERY AND/OR POST-PROCESSING) ! .............................................................. Use System_Constants ! For U_FLUX unit number IMPLICIT NONE INTEGER, INTENT(IN) :: ELEM_QP ! ELEM_QP = NUMBER OF QUADRATURE POINTS USED IN THIS LOOP IF ( U_FLUX > 0 ) WRITE (U_FLUX) ELEM_QP END SUBROUTINE SAVE_QP_LOOP_COUNT SUBROUTINE SAVE_QP_LOOP_ARRAYS (XYZ, E, B) ! .............................................................. ! SAVE ELEM POSITION, CONSTITUTIVE, & OPERATOR MATRICES AT QP ! USING (FOR SCP RECOVERY AND/OR POST-PROCESSING) ! .............................................................. Use System_Constants ! For U_FLUX, N_SPACE, N_R_B Use Elem_Type_Data ! for: LT_FREE IMPLICIT NONE REAL(DP), INTENT(IN) :: XYZ (N_SPACE) REAL(DP), INTENT(IN) :: E (N_R_B, N_R_B) REAL(DP), INTENT(IN) :: B (N_R_B, LT_FREE) IF ( U_FLUX > 0 ) WRITE (U_FLUX) XYZ, E, B END SUBROUTINE SAVE_QP_LOOP_ARRAYS SUBROUTINE GET_QP_LOOP_COUNT (ELEM_QP) ! .............................................................. ! GET NUMBER OF QUADRATURE POINTS THAT AN ELEMENT IS ! USING (FOR SCP RECOVERY AND/OR POST-PROCESSING) ! .............................................................. Use System_Constants ! For U_FLUX unit number IMPLICIT NONE INTEGER, INTENT(OUT) :: ELEM_QP ! ELEM_QP = NUMBER OF QUADRATURE POINTS USED IN THIS LOOP IF ( U_FLUX > 0 ) THEN READ (U_FLUX) ELEM_QP ELSE ELEM_QP = 0 PRINT *, 'GET_QP_LOOP_COUNT: WARNING, U_FLUX IS ZERO' N_WARN = N_WARN + 1 END IF END SUBROUTINE GET_QP_LOOP_COUNT SUBROUTINE GET_QP_LOOP_ARRAYS (XYZ, E, B) ! .............................................................. ! GET ELEM POSITION, CONSTITUTIVE, & OPERATOR MATRICES AT QP ! USING (FOR SCP RECOVERY AND/OR POST-PROCESSING) ! .............................................................. Use System_Constants ! For U_FLUX, N_SPACE, N_R_B Use Elem_Type_Data ! for: LT_FREE IMPLICIT NONE REAL(DP), INTENT(OUT) :: XYZ (N_SPACE) REAL(DP), INTENT(OUT) :: E (N_R_B, N_R_B) REAL(DP), INTENT(OUT) :: B (N_R_B, LT_FREE) IF ( U_FLUX > 0 ) THEN READ (U_FLUX) XYZ, E, B ELSE XYZ = 0.d0 ; B = 0.d0 CALL REAL_IDENTITY (N_R_B, E) PRINT *, 'GET_QP_LOOP_ARRAYS: WARNING, U_FLUX IS ZERO' N_WARN = N_WARN + 1 END IF END SUBROUTINE GET_QP_LOOP_ARRAYS SUBROUTINE GET_NODE_AT_DOF (N_DOF, N_G_FREE, N_PT) ! .............................................................. ! IDENTIFY NODE POINT, N_PT, ASSOCIATED WITH DEGREE OF FREEDOM ! .............................................................. IMPLICIT NONE INTEGER, INTENT(IN) :: N_DOF, N_G_FREE INTEGER, INTENT(OUT) :: N_PT ! N_DOF = EQUATION NUMBER ! N_G_FREE = NUMBER OF UNKNOWNS PER NODE ! N_PT = GLOBAL NODE NUMBER WHERE EQUATION OCCURS N_PT = 1 + (N_DOF -1)/N_G_FREE END SUBROUTINE GET_NODE_AT_DOF SUBROUTINE STORE_FLUX_POINT_COUNT ! Save LT_QP ! .............................................................. ! Store element flux point count ! .............................................................. Use System_Constants ! For U_FLUX, N_WARN Use Elem_Type_Data ! For LT_QP LOGICAL, SAVE :: WARN = .FALSE. IF ( U_FLUX > 0 ) THEN WRITE (U_FLUX) LT_QP ! Declare el flux save loop size ELSE IF ( NO_SCP_AVE ) RETURN IF ( .NOT. WARN ) THEN PRINT *,'WARNING: STORE_FLUX_POINT_COUNT, FILE NOT ACTIVE' PRINT *,'IS no_scp_ave IN KEYWORD CONTROLS ?' N_WARN = N_WARN + 1 WARN = .TRUE. END IF END IF END SUBROUTINE STORE_FLUX_POINT_COUNT ! Save LT_QP SUBROUTINE STORE_FLUX_POINT_DATA (XYZ, E, B) ! Save LT_QP ! .............................................................. ! Store element flux point data for averaging or post processing ! .............................................................. Use System_Constants ! For U_FLUX, N_WARN, N_R_B, N_SPACE Use Elem_Type_Data ! For LT_FREE REAL(DP), INTENT(IN) :: XYZ (N_SPACE), E (N_R_B, N_R_B), & B (N_R_B, LT_FREE) LOGICAL, SAVE :: WARN = .FALSE. IF ( U_FLUX > 0 ) THEN !--> WRITE COORDS, E AND DERIVATIVE MATRIX, FOR POST PROCESSING WRITE (U_FLUX) XYZ, E, B ELSE IF ( NO_SCP_AVE ) RETURN IF ( .NOT. WARN ) THEN PRINT *,'WARNING: STORE_FLUX_POINT_DATA, FILE NOT ACTIVE' PRINT *,'IS no_scp_ave IN KEYWORD CONTROLS ?' N_WARN = N_WARN + 1 WARN = .TRUE. END IF END IF END SUBROUTINE STORE_FLUX_POINT_DATA ! Save LT_QP SUBROUTINE READ_FLUX_POINT_COUNT (N_IP) ! Save LT_QP ! .............................................................. ! Read element flux point count ! .............................................................. Use System_Constants ! For U_FLUX, N_WARN Use Elem_Type_Data ! For LT_QP INTEGER, INTENT(OUT) :: N_IP LOGICAL, SAVE :: WARN = .FALSE. IF ( U_FLUX > 0 ) THEN READ (U_FLUX) N_IP ! Declare el flux save loop size ELSE N_IP = 0 IF ( NO_SCP_AVE ) RETURN IF ( .NOT. WARN ) THEN PRINT *,'WARNING: READ_FLUX_POINT_COUNT, FILE NOT ACTIVE' PRINT *,'IS no_scp_ave IN KEYWORD CONTROLS ?' N_WARN = N_WARN + 1 WARN = .TRUE. END IF END IF END SUBROUTINE READ_FLUX_POINT_COUNT ! Save LT_QP SUBROUTINE READ_FLUX_POINT_DATA (XYZ, E, B) ! Save LT_QP ! .............................................................. ! Read element flux point data for averaging or post processing ! .............................................................. Use System_Constants ! For U_FLUX, N_WARN, N_R_B, N_SPACE Use Elem_Type_Data ! For LT_FREE REAL(DP), INTENT(OUT) :: XYZ (N_SPACE), E (N_R_B, N_R_B), & B (N_R_B, LT_FREE) LOGICAL, SAVE :: WARN = .FALSE. IF ( U_FLUX > 0 ) THEN !--> READ COORDS, E AND DERIVATIVE MATRIX, FOR POST PROCESSING READ (U_FLUX) XYZ, E, B ELSE XYZ = 0.d0 ; E = 0.d0 ; B = 0.d0 IF ( NO_SCP_AVE ) RETURN IF ( .NOT. WARN ) THEN PRINT *,'WARNING: READ_FLUX_POINT_DATA, FILE NOT ACTIVE' PRINT *,'IS no_scp_ave IN KEYWORD CONTROLS ?' N_WARN = N_WARN + 1 WARN = .TRUE. END IF END IF END SUBROUTINE READ_FLUX_POINT_DATA ! Save LT_QP ! End File: model_lib.f