Merge pull request #139 from victorkemp/composite_buckling

Composite buckling

Author: zchlrnr

Source/EMG/EMG1/EMG.f90

@@ -97,6 +97,7 @@ SUBROUTINE EMG ( INT_ELEM_ID, OPT, WRITE_WARN, CALLING_SUBR, WRT_BUG_THIS_TIME )
 ! **********************************************************************************************************************************
       NUM_EMG_FATAL_ERRS = 0
 
+
       CALL EMG_INIT
 
 ! **********************************************************************************************************************************
@@ -343,7 +344,7 @@ SUBROUTINE EMG ( INT_ELEM_ID, OPT, WRITE_WARN, CALLING_SUBR, WRT_BUG_THIS_TIME )
          IF (NUM_EMG_FATAL_ERRS > 0)   CALL EMG_QUIT
 
       ELSE IF ((TYPE(1:5) == 'QUAD4') .OR. (TYPE == 'SHEAR   ')) THEN
-         CALL QDEL1 ( OPT, WRITE_WARN )
+         CALL QDEL1 ( OPT, INT_ELEM_ID, WRITE_WARN )
          IF (NUM_EMG_FATAL_ERRS > 0)   CALL EMG_QUIT
 
       ELSE IF ((TYPE == 'HEXA8   ') .OR. (TYPE == 'HEXA20  ') .OR.                                                                 &

Source/EMG/EMG4/QDEL1.f90

@@ -24,7 +24,7 @@
 
 ! End MIT license text.                                                                                      
 
-      SUBROUTINE QDEL1 ( OPT, WRITE_WARN )
+      SUBROUTINE QDEL1 ( OPT, INT_ELEM_ID, WRITE_WARN )
 
 ! Calculates, or calls subr's to calculate, quadrilateral element matrices:
 
@@ -54,6 +54,7 @@ SUBROUTINE QDEL1 ( OPT, WRITE_WARN )
 
       CHARACTER(LEN=LEN(BLNK_SUB_NAM)):: SUBR_NAME = 'QDEL1'
       CHARACTER(1*BYTE), INTENT(IN)   :: OPT(6)               ! 'Y'/'N' flags for whether to calc certain elem matrices
+      INTEGER(LONG), INTENT(IN)       :: INT_ELEM_ID          ! Internal element ID
       CHARACTER(LEN=*), INTENT(IN)    :: WRITE_WARN           ! If 'Y" write warning messages, otherwise do not
       CHARACTER( 1*BYTE)              :: RED_INT_SHEAR        ! If 'Y', use Gaussian weighted average of B matrices for shear terms
 
@@ -219,7 +220,7 @@ SUBROUTINE QDEL1 ( OPT, WRITE_WARN )
                ELSE
                   RED_INT_SHEAR = 'N'
                ENDIF
-               CALL QMEM1 ( OPT, IORD, RED_INT_SHEAR, AREA, XSD, YSD, BIG_BM )
+               CALL QMEM1 ( OPT, INT_ELEM_ID, IORD, RED_INT_SHEAR, AREA, XSD, YSD, BIG_BM )
             ENDIF
          ENDIF
 

Source/EMG/EMG4/QMEM1.f90

@@ -24,7 +24,7 @@
 
 ! End MIT license text.                                                                                      
 
-      SUBROUTINE QMEM1 ( OPT, IORD, RED_INT_SHEAR, AREA, XSD, YSD, BIG_BM )
+      SUBROUTINE QMEM1 ( OPT, INT_ELEM_ID, IORD, RED_INT_SHEAR, AREA, XSD, YSD, BIG_BM )
 
 ! Isoparametric membrane quadrilateral. Default iorq1s = 1 gives reduced integration for shear terms. User can override
 ! this with Bulk Data PARAM iorq1s 2. Element can be nonplanar. HBAR is the dist that the nodes are away from the mean
@@ -42,14 +42,16 @@ SUBROUTINE QMEM1 ( OPT, IORD, RED_INT_SHEAR, AREA, XSD, YSD, BIG_BM )
 
       USE PENTIUM_II_KIND, ONLY       :  BYTE, LONG, DOUBLE
       USE IOUNT1, ONLY                :  ERR, F04, F06, WRT_ERR, WRT_LOG
-      USE SCONTR, ONLY                :  BLNK_SUB_NAM, FATAL_ERR, MAX_ORDER_GAUSS, MAX_STRESS_POINTS, MEFE, NSUB, NTSUB
+      USE SCONTR, ONLY                :  BLNK_SUB_NAM, FATAL_ERR, MAX_ORDER_GAUSS, MAX_STRESS_POINTS, MEFE, NSUB, NTSUB,           &
+                                         MRPCOMP_PLIES, MRPCOMP0
       USE TIMDAT, ONLY                :  TSEC
       USE SUBR_BEGEND_LEVELS, ONLY    :  QMEM1_BEGEND
       USE CONSTANTS_1, ONLY           :  ZERO, FOUR
       USE NONLINEAR_PARAMS, ONLY      :  LOAD_ISTEP
       USE MODEL_STUF, ONLY            :  ALPVEC, BE1, BMEANT, DT, EID, ELDOF, ELGP, EM, ERR_SUB_NAM, HBAR, KE, KED, MXWARP,        &
-                                         NUM_EMG_FATAL_ERRS, PCOMP_LAM, PCOMP_PROPS, PPE, PRESS, PTE,                              &
-                                         SE1, STE1, SHELL_AALP, SHELL_A, TREF, TYPE, FCONV, STRESS
+                                         NUM_EMG_FATAL_ERRS, PCOMP_LAM, PCOMP_PROPS, RPCOMP, PPE, PRESS, PTE,                      &
+                                         SE1, STE1, SHELL_AALP, SHELL_A, TREF, TYPE, FCONV, STRESS, NUM_PLIES, INTL_PID, TPLY,     &
+                                         EPROP, PLY_NUM
       USE DEBUG_PARAMETERS, ONLY      :  DEBUG
 
       USE QMEM1_USE_IFs
@@ -61,9 +63,11 @@ SUBROUTINE QMEM1 ( OPT, IORD, RED_INT_SHEAR, AREA, XSD, YSD, BIG_BM )
       CHARACTER( 1*BYTE), INTENT(IN)  :: RED_INT_SHEAR     ! If 'Y', use Gaussian weighted average of B matrices for shear terms
       CHARACTER(46*BYTE)              :: IORD_MSG          ! Character name of the integration order (used for debug output)
 
+      INTEGER(LONG), INTENT(IN)       :: INT_ELEM_ID       ! Internal element ID
       INTEGER(LONG), INTENT(IN)       :: IORD              ! Gaussian integration order for element
       INTEGER(LONG)                   :: GAUSS_PT          ! Gauss point number (used for output in subr SHP2DQ
       INTEGER(LONG)                   :: I,J,K,L,M         ! DO loop indices
+      INTEGER(LONG)                   :: JPLY              ! PLY_NUM in a DO loop
       INTEGER(LONG), PARAMETER        :: ID1( 8) = (/ 1, & ! ID1(1) =  1 means virgin  8X8  elem DOF  1 is MYSTRAN 24X24 elem DOF  1
                                                       2, & ! ID1(2) =  2 means virgin  8X8  elem DOF  2 is MYSTRAN 24X24 elem DOF  2
                                                       7, & ! ID1(3) =  7 means virgin  8X8  elem DOF  3 is MYSTRAN 24X24 elem DOF  7
@@ -95,7 +99,8 @@ SUBROUTINE QMEM1 ( OPT, IORD, RED_INT_SHEAR, AREA, XSD, YSD, BIG_BM )
       INTEGER(LONG), PARAMETER        :: NUM_NODES = 4     ! Quad has 4 nodes
                                                            ! Indicator of no output of elem data to BUG file
       INTEGER(LONG), PARAMETER        :: SUBR_BEGEND = QMEM1_BEGEND
-  
+      INTEGER(LONG)                   :: PLY_RPCOMP_INDEX  ! Index in array RPCOMP where data for ply K begins
+
       REAL(DOUBLE) , INTENT(IN)       :: AREA              ! Element area
       REAL(DOUBLE) , INTENT(IN)       :: XSD(4)            ! Diffs in x coords of quad sides in local coords
       REAL(DOUBLE) , INTENT(IN)       :: YSD(4)            ! Diffs in y coords of quad sides in local coords
@@ -138,16 +143,16 @@ SUBROUTINE QMEM1 ( OPT, IORD, RED_INT_SHEAR, AREA, XSD, YSD, BIG_BM )
       REAL(DOUBLE)                    :: NBAR(2,8)         ! Matrix of shape functions (used in PPE calc)
       REAL(DOUBLE)                    :: PQ(8,NSUB)        ! Element pressure load matrix for the 8 local emenent DOF's.
       REAL(DOUBLE)                    :: PSH(4)            ! Output from subr SHP2DQ, called herein.
-      REAL(DOUBLE)                    :: FORCEx            ! Engineering force in the elem x direction
-      REAL(DOUBLE)                    :: FORCEy            ! Engineering force in the elem x direction
-      REAL(DOUBLE)                    :: FORCExy           ! Engineering force in the elem xy direction
+      REAL(DOUBLE)                    :: FORCEx(IORD_STRESS_Q4*IORD_STRESS_Q4) ! Engineering force in the elem x direction at Gauss points
+      REAL(DOUBLE)                    :: FORCEy(IORD_STRESS_Q4*IORD_STRESS_Q4) ! Engineering force in the elem x direction at Gauss points
+      REAL(DOUBLE)                    :: FORCExy(IORD_STRESS_Q4*IORD_STRESS_Q4)! Engineering force in the elem xy direction at Gauss points
       REAL(DOUBLE)                    :: QLOAD(2,NSUB)     ! 2 components of the element pressure load (from PRESS) in 1 array
       REAL(DOUBLE)                    :: SSS(MAX_ORDER_GAUSS)
                                                            ! An output from subr ORDER, called herein. Gauss abscissa's.
       REAL(DOUBLE)                    :: SUMB              ! An intermediate variable used in calc B matrix for reduced integration
       REAL(DOUBLE)                    :: SUMD              ! An intermediate variable used in calc B matrix for reduced integration
       REAL(DOUBLE)                    :: TBAR              ! Average elem temperature 
- 
+
 ! **********************************************************************************************************************************
       IF (WRT_LOG >= SUBR_BEGEND) THEN
          CALL OURTIM
@@ -525,6 +530,65 @@ SUBROUTINE QMEM1 ( OPT, IORD, RED_INT_SHEAR, AREA, XSD, YSD, BIG_BM )
 
       IF ((OPT(6) == 'Y') .AND. (LOAD_ISTEP > 1)) THEN
 
+        ! Things not tested or properly coded for differential stiffness matrix yet:
+        ! - Non-symmetric composite
+        ! - Offset
+        ! - Non-planar element
+        ! - Thermal loading with composite
+        ! - Differential stiffness solution type
+
+                                                           ! Find membrane engineering forces at each Gauss point
+                                                           ! by summing the forces in each ply.
+        DO GAUSS_PT=1,IORD_STRESS_Q4*IORD_STRESS_Q4
+          FORCEx(GAUSS_PT)  = 0
+          FORCEy(GAUSS_PT)  = 0
+          FORCExy(GAUSS_PT) = 0
+        ENDDO
+
+        CALL GET_ELEM_NUM_PLIES ( INT_ELEM_ID )            ! Get NUM_PLIES
+ 
+        DO JPLY=1,NUM_PLIES
+                                                           ! Get UEL, EM, TPLY for this ply.
+          IF (PCOMP_PROPS == 'N') THEN
+                                                           ! EM is already set above.
+
+            TPLY = EPROP(1)                                ! Element thickness
+            CALL ELMDIS
+
+          ELSE
+
+                                                           
+            IF (PCOMP_LAM == 'NON') THEN                   ! Delete this IF block to allow the LAM field set to nonsymmetric.
+              FATAL_ERR          = FATAL_ERR + 1
+              NUM_EMG_FATAL_ERRS = NUM_EMG_FATAL_ERRS + 1
+              WRITE(ERR,*) ' *ERROR: Code not written for non-symmetric composite buckling or differential stiffness.'
+              WRITE(F06,*) ' *ERROR: Code not written for non-symmetric composite buckling or differential stiffness.'
+              CALL OUTA_HERE ( 'Y' )
+            ENDIF
+
+            PLY_NUM = JPLY                                 ! Used by SHELL_ABD_MATRICES
+            CALL SHELL_ABD_MATRICES ( INT_ELEM_ID, 'N' )   ! Get EM, ZPLY, TPLY, ALPVEC for this ply
+            CALL ELMDIS
+            CALL ELMDIS_PLY                                ! Adjust UEL using ZPLY
+
+          ENDIF
+
+          GAUSS_PT = 0
+          DO I=1,IORD_STRESS_Q4
+            DO J=1,IORD_STRESS_Q4
+              GAUSS_PT = GAUSS_PT + 1
+
+              CALL ELEM_STRE_STRN_ARRAYS ( GAUSS_PT+1 )    ! Stress at this Gauss point using UEL, BE1, EM, ALPVEC, DT
+                                                             
+              FORCEx(GAUSS_PT)  = FORCEx(GAUSS_PT)  + TPLY*STRESS(1)
+              FORCEy(GAUSS_PT)  = FORCEy(GAUSS_PT)  + TPLY*STRESS(2)
+              FORCExy(GAUSS_PT) = FORCExy(GAUSS_PT) + TPLY*STRESS(3)
+            ENDDO
+          ENDDO
+          
+        ENDDO
+
+
 ! Accoring to:
 !   Robert D. Cook, David S. Malkus, Michael E. Plesha Concepts and Applications of Finite Element Analysis, 3rd Edition  1989
 !   Section 14.3 Stress Stiffness Matrix Of A Plate Element
@@ -547,28 +611,21 @@ SUBROUTINE QMEM1 ( OPT, IORD, RED_INT_SHEAR, AREA, XSD, YSD, BIG_BM )
 ! [ σ] = ⌡   ⌡   DPSHX  [ Nxy Ny  ]  DPSHX  |J|  dξ dη
 !        -1  -1  
 
-        CALL ELMDIS
-
-        CALL ORDER_GAUSS ( IORD_STRESS_Q4, SSS, HHH )
-
         DO K=1,ELGP
           DO L=1,ELGP
             KS(K,L) = ZERO
           ENDDO   
         ENDDO 
 
+        CALL ORDER_GAUSS ( IORD_STRESS_Q4, SSS, HHH )
+
         GAUSS_PT = 0
         DO I=1,IORD_STRESS_Q4
           DO J=1,IORD_STRESS_Q4
             GAUSS_PT = GAUSS_PT + 1
 
-            CALL ELEM_STRE_STRN_ARRAYS (GAUSS_PT+1)        ! Stress at this Gauss point
-                                                           
-            FORCEx  = FCONV(1)*STRESS(1)                   ! Engineering forces at this Gauss point
-            FORCEy  = FCONV(1)*STRESS(2)
-            FORCExy = FCONV(1)*STRESS(3)
-            DUM11(1,1) = FORCEx  ; DUM11(1,2) = FORCExy
-            DUM11(2,1) = FORCExy ; DUM11(2,2) = FORCEy
+            DUM11(1,1) = FORCEx(GAUSS_PT)  ; DUM11(1,2) = FORCExy(GAUSS_PT)
+            DUM11(2,1) = FORCExy(GAUSS_PT) ; DUM11(2,2) = FORCEy(GAUSS_PT)
 
             DO L=1,ELGP                                    ! Shape function derivatives at this Gauss point.
               DPSHX(1,L) = BE1(1,ID1(2*(L-1)+1),GAUSS_PT+1)     
@@ -589,6 +646,7 @@ SUBROUTINE QMEM1 ( OPT, IORD, RED_INT_SHEAR, AREA, XSD, YSD, BIG_BM )
 
           ENDDO
         ENDDO
+
                                                            ! Copy KS into KED for each translational DOF.
         DO I=1,6*ELGP
           DO J=1,6*ELGP
@@ -603,7 +661,6 @@ SUBROUTINE QMEM1 ( OPT, IORD, RED_INT_SHEAR, AREA, XSD, YSD, BIG_BM )
           ENDDO   
         ENDDO 
 
-
       ENDIF
 
 ! **********************************************************************************************************************************

Source/EMG/EMG4/QPLT2.f90

@@ -265,6 +265,7 @@ SUBROUTINE QPLT2 ( OPT, AREA, XSD, YSD, BIG_BB )
 
          GAUSS_PT = 0
          IORD_MSG = 'for plate bending strains,           IORQ2B = '
+
          DO I=1,IORQ2B
             DO J=1,IORQ2B
                SSI = SSS(I)
@@ -307,7 +308,6 @@ SUBROUTINE QPLT2 ( OPT, AREA, XSD, YSD, BIG_BB )
          ENDDO 
 
          IF (IORQ2T > 0) THEN
-
             GAUSS_PT = 0
             IORD_MSG = 'for transverse shear strains,        IORQ2T = '
             DO I=1,IORDXX
@@ -564,7 +564,7 @@ SUBROUTINE QPLT2 ( OPT, AREA, XSD, YSD, BIG_BB )
       ENDDO
 
       IF ((PCOMP_PROPS == 'Y') .AND. (PCOMP_LAM == 'NON')) THEN
-
+         CALL ORDER_GAUSS ( IORQ2B, SSS, HHH )
          GAUSS_PT = 0
          IORD_MSG = 'for plate bending strains,           IORQ2B = '
          DO I=1,IORQ2B

Source/EMG/EMG4/TMEM1.f90

@@ -37,14 +37,15 @@ SUBROUTINE TMEM1 ( OPT, AREA, X2E, X3E, Y3E, WRT_BUG_THIS_TIME, BIG_BM )
 !  5) KED       = element differen stiff matrix calc   , if OPT(6) = 'Y' = 'Y'
 
       USE PENTIUM_II_KIND, ONLY       :  BYTE, LONG, DOUBLE
-      USE IOUNT1, ONLY                :  BUG, F04, WRT_BUG, WRT_LOG, F06
-      USE SCONTR, ONLY                :  BLNK_SUB_NAM, ELDT_BUG_BCHK_BIT, ELDT_BUG_BMAT_BIT, NSUB, NTSUB
+      USE IOUNT1, ONLY                :  ERR, BUG, F04, WRT_BUG, WRT_LOG, F06
+      USE SCONTR, ONLY                :  BLNK_SUB_NAM, ELDT_BUG_BCHK_BIT, ELDT_BUG_BMAT_BIT, NSUB, NTSUB, FATAL_ERR
       USE TIMDAT, ONLY                :  TSEC
       USE SUBR_BEGEND_LEVELS, ONLY    :  TMEM1_BEGEND
       USE CONSTANTS_1, ONLY           :  ZERO, ONE, THREE
       USE NONLINEAR_PARAMS, ONLY      :  LOAD_ISTEP
       USE MODEL_STUF, ONLY            :  ALPVEC, BE1, EID, DT, EM, ELDOF, KE, PCOMP_LAM, PCOMP_PROPS, PRESS, PPE, PTE, SE1, STE1,  &
-                                         SHELL_AALP, SHELL_A, SHELL_PROP_ALP, TREF, TYPE, XEB, XEL, ELGP, FCONV, STRESS, KED
+                                         SHELL_AALP, SHELL_A, SHELL_PROP_ALP, TREF, TYPE, XEB, XEL, ELGP, FCONV, STRESS, KED,      &
+                                         NUM_EMG_FATAL_ERRS
       USE DEBUG_PARAMETERS, ONLY      :  DEBUG
 
       USE TMEM1_USE_IFs
@@ -262,6 +263,14 @@ SUBROUTINE TMEM1 ( OPT, AREA, X2E, X3E, Y3E, WRT_BUG_THIS_TIME, BIG_BM )
 
       IF ((OPT(6) == 'Y') .AND. (LOAD_ISTEP > 1)) THEN
 
+        IF (PCOMP_PROPS == 'Y') THEN
+          FATAL_ERR          = FATAL_ERR + 1
+          NUM_EMG_FATAL_ERRS = NUM_EMG_FATAL_ERRS + 1
+          WRITE(ERR,*) ' *ERROR: Code not written for CTRIA3 composite buckling or differential stiffness.'
+          WRITE(F06,*) ' *ERROR: Code not written for CTRIA3 composite buckling or differential stiffness.'
+          CALL OUTA_HERE ( 'Y' )
+        ENDIF
+
 ! Accoring to:
 !   Robert D. Cook, David S. Malkus, Michael E. Plesha Concepts and Applications of Finite Element Analysis, 3rd Edition  1989
 !   Section 14.3 Stress Stiffness Matrix Of A Plate Element

Source/Interfaces/QDEL1_Interface.f90

@@ -28,7 +28,7 @@ MODULE QDEL1_Interface
 
    INTERFACE
 
-      SUBROUTINE QDEL1 ( OPT, WRITE_WARN )
+      SUBROUTINE QDEL1 ( OPT, INT_ELEM_ID, WRITE_WARN )
 
 
       USE PENTIUM_II_KIND, ONLY       :  BYTE, LONG, DOUBLE
@@ -46,6 +46,7 @@ SUBROUTINE QDEL1 ( OPT, WRITE_WARN )
       IMPLICIT NONE 
 
       CHARACTER(1*BYTE), INTENT(IN)   :: OPT(6)               ! 'Y'/'N' flags for whether to calc certain elem matrices
+      INTEGER(LONG), INTENT(IN)       :: INT_ELEM_ID          ! Internal element ID
       CHARACTER(LEN=*), INTENT(IN)    :: WRITE_WARN           ! If 'Y" write warning messages, otherwise do not
       CHARACTER( 1*BYTE)              :: RED_INT_SHEAR        ! If 'Y', use Gaussian weighted average of B matrices for shear terms
 

Source/Interfaces/QMEM1_Interface.f90

@@ -28,7 +28,7 @@ MODULE QMEM1_Interface
 
    INTERFACE
 
-      SUBROUTINE QMEM1 ( OPT, IORD, RED_INT_SHEAR, AREA, XSD, YSD, BIG_BM )
+      SUBROUTINE QMEM1 ( OPT, INT_ELEM_ID, IORD, RED_INT_SHEAR, AREA, XSD, YSD, BIG_BM )
 
 
       USE PENTIUM_II_KIND, ONLY       :  BYTE, LONG, DOUBLE
@@ -48,6 +48,7 @@ SUBROUTINE QMEM1 ( OPT, IORD, RED_INT_SHEAR, AREA, XSD, YSD, BIG_BM )
       CHARACTER( 1*BYTE), INTENT(IN)  :: RED_INT_SHEAR     ! If 'Y', use Gaussian weighted average of B matrices for shear terms
       CHARACTER(46*BYTE)              :: IORD_MSG          ! Character name of the integration order (used for debug output)
 
+      INTEGER(LONG), INTENT(IN)       :: INT_ELEM_ID       ! Internal element ID
       INTEGER(LONG), INTENT(IN)       :: IORD              ! Gaussian integration order for element
       INTEGER(LONG)                   :: GAUSS_PT          ! Gauss point number (used for output in subr SHP2DQ
       INTEGER(LONG), PARAMETER        :: ID1( 8) = (/ 1, & ! ID1(1) =  1 means virgin  8X8  elem DOF  1 is MYSTRAN 24X24 elem DOF  1

Source/LK1/L1A/LOADB.f90

@@ -615,8 +615,8 @@ SUBROUTINE LOADB
       !  If they want to continue put DEBUG 201 w/val /= 0 entry in the BDF
 
       IF ((SOL_NAME(1:8) == 'DIFFEREN') .OR. (SOL_NAME(1:8) == 'BUCKLING')) THEN
-         IF((NCTRIA3   > 0) .OR. (NCTRIA3K  > 0) .OR. (NCQUAD4   > 0) .OR. (NCQUAD4K  > 0) .OR. (NCSHEAR   > 0)) THEN
-            MESSAG1= ' *WARNING: BUCKLING and DIFFERN SOL are only coded for the BAR, HEXA, PENTA, and TETRA elements'
+         IF((NCSHEAR   > 0)) THEN
+            MESSAG1= ' *WARNING: BUCKLING and DIFFERN SOL are only coded for the BAR, HEXA, PENTA, TETRA, TRIA, and QUAD elements'
             MESSAG2= '           Either remove all other elements or include a Bulk Data entry: DEBUG   201, with value /= 0'
             IF (DEBUG(201) == 0) THEN 
                WRITE(F06,*) MESSAG1
@@ -642,8 +642,8 @@ SUBROUTINE LOADB
              ' ++',127X,'++',/,                                                                                                    &
    ' ++', 6X,' B U C K L I N G   &   D I F F E R E N   S O L   a r e   o n l y   c o d e d   f o r   e l e m e n t s:',18X,'++',/  &
              ' ++',127X,'++',/,                                                                                                    &
-   ' ++', 6X,'                                            BAR,  HEXA,  PENTA,  TETRA', 51X, '++'                                   &
-   8X,'++',/,' ++',127X,'++',/                                                                                                     &
+             ' ++                                BAR,  HEXA,  PENTA,  TETRA,  TRIA3,  TRIA3K,  QUAD4,  QUAD4K       ' ,28X,'++',/  &
+             ' ++',127X,'++',/                                                                                                     &
    ' ++', 6X,'   A l l   o t h e r   e l e m s   i g n o r e d   i n   t h e  d i f f e r e n t i a l   s t i f f   c a l c s',    &
   10X,'++',/,' ++',127X,'++',/                                                                                                     &
              ' ++',127X,'++',/                                                                                                     &