diff --git a/Source/EMG/EMG3/ELAS1.f90 b/Source/EMG/EMG3/ELAS1.f90 index 89998a76..a9c92aa5 100644 --- a/Source/EMG/EMG3/ELAS1.f90 +++ b/Source/EMG/EMG3/ELAS1.f90 @@ -75,11 +75,11 @@ SUBROUTINE ELAS1 ( OPT, WRITE_WARN ) ENDIF ! ********************************************************************************************************************************** -! Calculate SE1 matrix for stress recovery. +! Calculate SE1 matrix for force recovery. IF (OPT(3) == 'Y') THEN - SE1(1,I1,1) = K*FCONV(1) - SE1(1,I2,1) = -K*FCONV(1) + SE1(1,I1,1) = K + SE1(1,I2,1) = -K ENDIF diff --git a/Source/LK9/L92/ELEM_STRE_STRN_ARRAYS.f90 b/Source/LK9/L92/ELEM_STRE_STRN_ARRAYS.f90 index 263ae5b6..aa1b851a 100644 --- a/Source/LK9/L92/ELEM_STRE_STRN_ARRAYS.f90 +++ b/Source/LK9/L92/ELEM_STRE_STRN_ARRAYS.f90 @@ -98,6 +98,8 @@ SUBROUTINE ELEM_STRE_STRN_ARRAYS ( STR_PT_NUM ) ! ********************************************************************************************************************************** ! Calc stresses for 1D elements +! Warning: For ELAS1/2/3/4, the STRESS() array contains force, not stress. This is because we must calculate stress from force, not +! the other way around. IF ((TYPE(1:3) == 'BAR') .OR. (TYPE(1:4) == 'BUSH') .OR. (TYPE(1:4) == 'ELAS') .OR. (TYPE(1:3) == 'ROD') .OR. & (TYPE(1:5) == 'USER1')) THEN diff --git a/Source/LK9/L92/OFP3_ELFE_1D.f90 b/Source/LK9/L92/OFP3_ELFE_1D.f90 index bed546df..ffac7a5a 100644 --- a/Source/LK9/L92/OFP3_ELFE_1D.f90 +++ b/Source/LK9/L92/OFP3_ELFE_1D.f90 @@ -38,7 +38,7 @@ SUBROUTINE OFP3_ELFE_1D ( JVEC, FEMAP_SET_ID, ITE, OT4_EROW ) USE FEMAP_ARRAYS, ONLY : FEMAP_EL_NUMS, FEMAP_EL_VECS USE PARAMS, ONLY : OTMSKIP, PRTNEU USE MODEL_STUF, ONLY : ANY_ELFE_OUTPUT, BUSH_CID, BUSH_VVEC, EDAT, ELAS_COMP, ELEM_LEN_12, ELEM_LEN_AB, EPNT, & - ETYPE, EID, ELMTYP, ELOUT, FCONV, METYPE, NUM_EMG_FATAL_ERRS, OFFDIS_GA_GB, OFFDIS_L, & + ETYPE, EID, ELMTYP, ELOUT, METYPE, NUM_EMG_FATAL_ERRS, OFFDIS_GA_GB, OFFDIS_L, & PE_GA_GB, PEL, PLY_NUM, STRESS, TE, TE_GA_GB, TYPE, XEL USE LINK9_STUFF, ONLY : EID_OUT_ARRAY, MAXREQ, OGEL USE OUTPUT4_MATRICES, ONLY : OTM_ELFE, TXT_ELFE @@ -188,11 +188,7 @@ SUBROUTINE OFP3_ELFE_1D ( JVEC, FEMAP_SET_ID, ITE, OT4_EROW ) CALL ELEM_STRE_STRN_ARRAYS ( 1 ) NDUM = 0 CALL CALC_ELEM_STRESSES ( 1, NDUM, 0, 'N', 'N' ) - IF (FCONV(1) > ZERO) THEN ! ELAS engr force is stress/FCONV - OGEL(NUM_OGEL,1) = STRESS(1)/FCONV(1) - ELSE - OGEL(NUM_OGEL,1) = ZERO - ENDIF + OGEL(NUM_OGEL,1) = STRESS(1) ! ELAS engr force is stored in the stress array ! end elas ! --------------------------------------------------------------------------------------------------------------- ELSE IF (ETYPE(J)(1:4) == 'BUSH') THEN @@ -553,8 +549,6 @@ SUBROUTINE OFP3_ELFE_1D ( JVEC, FEMAP_SET_ID, ITE, OT4_EROW ) ENDIF CALL DEALLOCATE_FEMAP_DATA -! For ELAS we need to calculate elem engr forces from the stresses since there is no "local" elem coord system - ! elas1 --------------------------------------------------------------------------------------------------------------------------- NDUM = 0 NUM_FROWS= 0 @@ -577,11 +571,7 @@ SUBROUTINE OFP3_ELFE_1D ( JVEC, FEMAP_SET_ID, ITE, OT4_EROW ) CALL ELMDIS CALL ELEM_STRE_STRN_ARRAYS ( 1 ) CALL CALC_ELEM_STRESSES ( NCELAS1, NDUM, NUM_FROWS, 'N', 'Y' ) - IF (FCONV(1) > 0.D0) THEN - FEMAP_EL_VECS(NUM_FROWS,1) = STRESS(1)/FCONV(1) - ELSE - - ENDIF + FEMAP_EL_VECS(NUM_FROWS,1) = STRESS(1) ENDIF ENDDO IF (NUM_FROWS > 0) THEN @@ -611,11 +601,7 @@ SUBROUTINE OFP3_ELFE_1D ( JVEC, FEMAP_SET_ID, ITE, OT4_EROW ) CALL ELMDIS CALL ELEM_STRE_STRN_ARRAYS ( 1 ) CALL CALC_ELEM_STRESSES ( NCELAS2, NDUM, NUM_FROWS, 'N', 'Y' ) - IF (FCONV(1) > 0.D0) THEN - FEMAP_EL_VECS(NUM_FROWS,1) = STRESS(1)/FCONV(1) - ELSE - - ENDIF + FEMAP_EL_VECS(NUM_FROWS,1) = STRESS(1) ENDIF ENDDO IF (NUM_FROWS > 0) THEN @@ -645,11 +631,7 @@ SUBROUTINE OFP3_ELFE_1D ( JVEC, FEMAP_SET_ID, ITE, OT4_EROW ) CALL ELMDIS CALL ELEM_STRE_STRN_ARRAYS ( 1 ) CALL CALC_ELEM_STRESSES ( NCELAS3, NDUM, NUM_FROWS, 'N', 'Y' ) - IF (FCONV(1) > 0.D0) THEN - FEMAP_EL_VECS(NUM_FROWS,1) = STRESS(1)/FCONV(1) - ELSE - - ENDIF + FEMAP_EL_VECS(NUM_FROWS,1) = STRESS(1) ENDIF ENDDO IF (NUM_FROWS > 0) THEN @@ -679,11 +661,7 @@ SUBROUTINE OFP3_ELFE_1D ( JVEC, FEMAP_SET_ID, ITE, OT4_EROW ) CALL ELMDIS CALL ELEM_STRE_STRN_ARRAYS ( 1 ) CALL CALC_ELEM_STRESSES ( NCELAS4, NDUM, NUM_FROWS, 'N', 'Y' ) - IF (FCONV(1) > 0.D0) THEN - FEMAP_EL_VECS(NUM_FROWS,1) = STRESS(1)/FCONV(1) - ELSE - - ENDIF + FEMAP_EL_VECS(NUM_FROWS,1) = STRESS(1) ENDIF ENDDO IF (NUM_FROWS > 0) THEN diff --git a/Source/LK9/L92/ONE_D_STRESS_OUTPUTS.f90 b/Source/LK9/L92/ONE_D_STRESS_OUTPUTS.f90 index 9ba618a5..53487e4c 100644 --- a/Source/LK9/L92/ONE_D_STRESS_OUTPUTS.f90 +++ b/Source/LK9/L92/ONE_D_STRESS_OUTPUTS.f90 @@ -33,9 +33,8 @@ SUBROUTINE ONE_D_STRESS_OUTPUTS ( SIZE_ALLOCATED, NUM1, NUM_FEMAP_ROWS, WRITE_OG USE PENTIUM_II_KIND, ONLY : BYTE, LONG, DOUBLE USE IOUNT1, ONLY : ERR, F06 USE SCONTR, ONLY : BLNK_SUB_NAM, FATAL_ERR - USE TIMDAT, ONLY : TSEC USE CONSTANTS_1, ONLY : ZERO - USE MODEL_STUF, ONLY : STRESS, TYPE, ZS + USE MODEL_STUF, ONLY : STRESS, TYPE, ZS, FCONV USE LINK9_STUFF, ONLY : MSPRNT, OGEL USE FEMAP_ARRAYS, ONLY : FEMAP_EL_VECS USE PARAMS, ONLY : PRTNEU @@ -184,11 +183,11 @@ SUBROUTINE ONE_D_STRESS_OUTPUTS ( SIZE_ALLOCATED, NUM1, NUM_FEMAP_ROWS, WRITE_OG FATAL_ERR = FATAL_ERR + 1 CALL OUTA_HERE ( 'Y' ) ENDIF - OGEL(NUM1,1) = STRESS(1) + OGEL(NUM1,1) = FCONV(1)*STRESS(1) ! STRESS() contains force and FCONV(1) is S. ENDIF IF (WRITE_NEU .AND. (WRITE_FEMAP == 'Y')) THEN - FEMAP_EL_VECS(NUM_FEMAP_ROWS,1) = STRESS(1) - FEMAP_EL_VECS(NUM_FEMAP_ROWS,2) = STRESS(1) + FEMAP_EL_VECS(NUM_FEMAP_ROWS,1) = FCONV(1)*STRESS(1) + FEMAP_EL_VECS(NUM_FEMAP_ROWS,2) = FCONV(1)*STRESS(1) ENDIF ELSE IF (TYPE == 'ROD ') THEN ! ROD1 elements diff --git a/Source/Modules/MODEL_STUF.f03 b/Source/Modules/MODEL_STUF.f03 index 363116b5..13dc279d 100644 --- a/Source/Modules/MODEL_STUF.f03 +++ b/Source/Modules/MODEL_STUF.f03 @@ -1421,6 +1421,7 @@ MODULE MODEL_STUF REAL(DOUBLE) :: FCONV(3) = (/(ZERO, I=1,3)/) ! Array of constants to convert stresses to engr forces for elem + ! Except CELAS where it converts force to stress REAL(DOUBLE) :: HBAR = ZERO ! For the quad elem, the dist from the mean plane to the G.P.'s @@ -1511,6 +1512,7 @@ MODULE MODEL_STUF REAL(DOUBLE) :: STRESS(9) = (/(ZERO, I=1,9)/) ! Array of elem stresses for one S/C + ! Except CELAS where it is force, not stress. REAL(DOUBLE) :: TE(3,3) = RESHAPE ( (/(ZERO, I=1,3*3)/), (/3,3/) ) ! Coord system transformation matrix (3 x 3) such that UEL = TE*UEB