! copyright 2005,  J. E. Akin, all rights reserved.
! Begin application_lib.f

! NOTE: This is the only set of programs that vary from one application
! to the next.  Only ELEM_SQ_MATRIX is always required by every problem.
! It must be edited to contain the problem dependent source, or the 
! include path to the problem dependent source.  Then this library must
! be compiled and linked to the other archived executables.

SUBROUTINE 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)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * *
! GENERATE ELEMENT SQUARE MATRIX, OPTIONAL COLUMN MATRIX
! * * * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants
Use Geometric_Properties 
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_M (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 Sys_Properties_Data  ! for FLT_MISC, MISC_FIX
Use Interface_Header     ! for functions
Use Select_source        ! for analytic data
 IMPLICIT NONE

 INTEGER,  INTENT(IN)    :: IE
 REAL(DP), INTENT(INOUT) :: E (N_R_B, N_R_B)
 REAL(DP), INTENT(INOUT) :: H_INTG (LT_N)     !b make optional !uu
 REAL(DP), INTENT(IN)    :: X (MAX_NP, N_SPACE) ! Global coordinates

!    OPTIONAL PROPERTY VALUES
 REAL(DP), INTENT(IN) :: PRT_L_PT (LT_N, N_NP_FLO),              &
                         PRT_MAT (MISC_FL) 
 INTEGER,  INTENT(IN) :: L_PT_PROP (LT_N, N_NP_FIX)
!b XXX Error L_PT_PROP (LT_N, N_NP_FIX) is correct XXX
                                                                   
!         Automatic Arrays
 REAL(DP) :: AJ (N_SPACE, N_SPACE),     XYZ (N_SPACE),           &
             AJ_INV (N_SPACE, N_SPACE), BODY (N_SPACE)
 REAL(DP) :: STRAIN_0 (N_R_B),          STRAIN (N_R_B + 2),      &
             STRESS (N_R_B + 2)
 !b REAL(DP) :: H_INTG (LT_N)
 REAL(DP) :: B (N_R_B, LT_FREE), EB (N_R_B, LT_FREE),            &
             DGH (N_SPACE, LT_N), DGV (N_SPACE, LT_FREE)

!                 VARIABLES:
! AJ        = JACOBIAN
! AJ_INV    = JACOBIAN INVERSE
! B         = STRAIN-DISPLACEMENT (GRADIENT) MATRIX
! BODY      = BODY FORCE VECTOR
! COORD     = SPATIAL COORDINATES OF ELEMENT'S NODES
! D         = NODAL PARAMETERS ASSOCIATED WITH AN ELEMENT
! DGH       = GLOBAL DERIVATIVES SCALAR INTERPOLATION FUNCTIONS
! DGV       = GLOBAL DERIVATIVES VECTOR INTERPOLATION FUNCTIONS
! DIAG_M    = ELEMENT DIAGONAL MASS MATRIX
! DLG       = LOCAL DERIVATIVES GEOMETRIC INTERPOLATION
! DLH       = LOCAL DERIVATIVES SCALAR INTERPOLATION FUNCTIONS
! DLV       = LOCAL DERIVATIVES VECTOR INTERPOLATION FUNCTIONS
! E         = CONSTITUTIVE MATRIX
! EB        = PRODUCT OF E*B
! EL_M      = ELEMENT CONSISTENT MASS MATRIX
! FLT_MISC  = SYSTEM STORAGE OF FLOATING PT MISC PROP
! G         = GEOMETRIC INTERPOLATION FUNCTIONS
! H         = SCALAR INTERPOLATION FUNCTIONS
! H_INTG    = INTEGRAL OF INTERPOLATION FUNCTIONS
! L_PT_PROP = INTEGER PROPERTIES AT EACH ELEMENT NODE
! MISC_FIX  = MISCELLANEOUS INTEGER SYSTEM PROPERTIES
! LT_N      = NUMBER OF NODES PER ELEMENT
! LT_FREE   = NUMBER OF DEGREES OF FREEDOM PER ELEMENT
! LT_GEOM   = NUMBER OF GEOMETRY NODES
! N_MAT     = NUMBER OF MATERIAL TYPES
! LT_PARM   = DIMENSION OF PARAMETRIC SPACE
! 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   = UNIT FOR POST SOLUTION MATRICES STORAGE
! PRT_L_PT  = REAL PROPERTIES AT ELEMENT NODES
! PRT_MAT   = REAL ELEM PROPERTIES BASED ON MATERIAL NUMBER
! PT        = QUADRATURE COORDINATES
! S         = ELEMENT SQUARE MATRIX
! STRAIN    = STRAIN OR GRADIENT VECTOR (PLUS TWO EXTRA)
! STRAIN_0  = INITIAL STRAIN OR GRADIENT VECTOR
! STRESS    = STRESS VECTOR (PLUS TWO EXTRA)
! V         = VECTOR INTERPOLATION FUNCTIONS
! WT        = QUADRATURE WEIGHTS
! XYZ       = SPACE COORDINATES AT A POINT
! ..............................................................
!    ***  ELEM_SQ_MATRIX PROBLEM DEPENDENT STATEMENTS FOLLOW ***
! ..............................................................
  include 'my_el_sq_inc' ! MUST use, or hardcode S here, C optional

!     Allow a CONTAINS statement (no statements after include)
!  IF ( DEBUG_EL_SQ ) PRINT *, 'Leaving EL_SQ ', THIS_EL
! suppress compiler warnings (touch is never true)
 !IF ( TOUCH ) CALL TOUCH_UNUSED_ITEMS_1 (AJ, AJ_INV, B, BODY, DGH, &
 !             DGV, E, EB, G, H_INTG, IE, L_PT_PROP, PRT_L_PT,      &
 !             PRT_MAT, STRAIN, STRAIN_0, STRESS, XYZ)

END SUBROUTINE ELEM_SQ_MATRIX
                                   
SUBROUTINE MIXED_SQ_MATRIX (E, H_INTG, PRT_L_PT, PRT_MAT, &
                            L_PT_PROP, IE)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * *
! GENERATE MIXED_BC SQUARE MATRIX, OPTIONAL COLUMN MATRIX
! * * * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants
Use Geometric_Properties 
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 Sys_Properties_Data  ! for FLT_MISC, MISC_FIX
Use Interface_Header     ! for functions
Use Select_Source
 IMPLICIT NONE

 INTEGER,  INTENT(IN)    :: IE

!    OPTIONAL PROPERTY VALUES
 REAL(DP), INTENT(INOUT) :: E (N_R_B, N_R_B)
 REAL(DP), INTENT(IN)    :: PRT_L_PT (LT_N, N_NP_FLO),           &
                            PRT_MAT (MISC_FL) 
 INTEGER,  INTENT(IN)    :: L_PT_PROP (LT_N, N_NP_FIX)
                                                                   
!         Automatic Arrays
 REAL(DP) :: AJ (N_SPACE, N_SPACE),     XYZ (N_SPACE),           &
             AJ_INV (N_SPACE, N_SPACE), BODY (N_SPACE)
 REAL(DP) :: STRAIN_0 (N_R_B),          STRAIN (N_R_B + 2),      &
             STRESS (N_R_B + 2)
 REAL(DP) :: B (N_R_B, LT_FREE), EB (N_R_B, LT_FREE), H_INTG (LT_N), &
             DGH (N_SPACE, LT_N), DGV (N_SPACE, LT_FREE)

!                 VARIABLES:
! AJ        = JACOBIAN
! AJ_INV    = JACOBIAN INVERSE
! B         = STRAIN-DISPLACEMENT (GRADIENT) MATRIX
! BODY      = BODY FORCE VECTOR
! COORD     = SPATIAL COORDINATES OF ELEMENT'S NODES
! D         = NODAL PARAMETERS ASSOCIATED WITH AN ELEMENT
! 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
! E         = CONSTITUTIVE MATRIX
! EB        = PRODUCT OF E*B
! FLT_MISC  = SYSTEM STORAGE OF FLOATING PT MISC PROP
! G         = GEOMETRIC INTERPOLATION FUNCTIONS
! H         = SCALAR INTERPOLATION FUNCTIONS
! H_INTG    = INTEGRAL OF INTERPOLATION FUNCTIONS
! L_PT_PROP = INTEGER PROPERTIES AT EACH ELEMENT NODE
! MISC_FIX  = MISCELLANEOUS INTEGER SYSTEM PROPERTIES
! LT_N      = NUMBER OF NODES PER ELEMENT
! LT_FREE   = NUMBER OF DEGREES OF FREEDOM PER ELEMENT
! LT_GEOM   = NUMBER OF GEOMETRY NODES
! N_MAT     = NUMBER OF MATERIAL TYPES
! LT_PARM   = DIMENSION OF PARAMETRIC SPACE
! 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   = UNIT FOR POST SOLUTION MATRICES STORAGE
! PRT_L_PT  = REAL PROPERTIES AT ELEMENT NODES
! PRT_MAT   = REAL ELEM PROPERTIES BASED ON MATERIAL NUMBER
! PT        = QUADRATURE COORDINATES
! S         = ELEMENT SQUARE MATRIX
! STRAIN    = STRAIN OR GRADIENT VECTOR (PLUS TWO EXTRA)
! STRAIN_0  = INITIAL STRAIN OR GRADIENT VECTOR
! STRESS    = STRESS VECTOR (PLUS TWO EXTRA)
! V         = VECTOR INTERPOLATION FUNCTIONS
! WT        = QUADRATURE WEIGHTS
! XYZ       = SPACE COORDINATES AT A POINT
! ..............................................................
!    ***  MIXED_SQ_MATRIX PROBLEM DEPENDENT STATEMENTS FOLLOW ***
! ..............................................................

  include 'my_mixed_sq_inc' ! MUST use, or hardcode S here, C optional

!     Allow a CONTAINS statement (no statements after include)
! suppress compiler warnings (never true)
! IF ( TOUCH ) CALL TOUCH_UNUSED_ITEMS_1 (AJ, AJ_INV, B, BODY, DGH, &
!              DGV, E, EB, G, H_INTG, IE, L_PT_PROP, PRT_L_PT,      &
!              PRT_MAT, STRAIN, STRAIN_0, STRESS, XYZ) 
END SUBROUTINE MIXED_SQ_MATRIX
                                   
SUBROUTINE ELEM_COL_MATRIX (E, H_INTG, PRT_L_PT, PRT_MAT,&
                            L_PT_PROP, IE)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * *
!        GENERATE OPTIONAL ELEMENT COLUMN MATRIX
! * * * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants
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 Sys_Properties_Data  ! for FLT_MISC, MISC_FIX
Use Interface_Header     ! for functions
Use Select_source        ! for analytic data
 IMPLICIT NONE

 INTEGER,  INTENT(IN) :: IE

!    OPTIONAL PROPERTY VALUES
 REAL(DP), INTENT(INOUT) :: E (N_R_B, N_R_B)
 REAL(DP), INTENT(IN)    :: PRT_L_PT (LT_N, N_NP_FLO), &
                            PRT_MAT (MISC_FL) 
 INTEGER,  INTENT(IN)    :: L_PT_PROP (LT_N, N_NP_FIX)

!         Automatic Arrays
 REAL(DP) :: AJ (N_SPACE, N_SPACE), AJ_INV (N_SPACE, N_SPACE), &
             XYZ (N_SPACE), BODY (N_SPACE)
 REAL(DP) :: STRAIN_0 (N_R_B),          STRAIN (N_R_B + 2),      &
             STRESS (N_R_B + 2)
 REAL(DP) :: B (N_R_B, LT_FREE), EB (N_R_B, LT_FREE), H_INTG (LT_N), &
             DGH (N_SPACE, LT_N), DGV (N_SPACE, LT_FREE)
 
!                 VARIABLES:
! AJ        = JACOBIAN
! AJ_INV    = JACOBIAN INVERSE
! B         = STRAIN-DISPLACEMENT (GRADIENT) MATRIX
! BODY      = BODY FORCE VECTOR
! COORD     = SPATIAL COORDINATES OF ELEMENT'S NODES
! D         = NODAL PARAMETERS ASSOCIATED WITH AN ELEMENT
! 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
! E         = CONSTITUTIVE MATRIX
! EB        = PRODUCT OF E*B
! FLT_MISC  = SYSTEM STORAGE OF FLOATING PT MISC PROP
! G         = GEOMETRIC INTERPOLATION FUNCTIONS
! H         = SCALAR INTERPOLATION FUNCTIONS
! H_INTG    = INTEGRAL OF INTERPOLATION FUNCTIONS
! L_PT_PROP = INTEGER PROPERTIES AT EACH ELEMENT NODE
! MISC_FIX  = MISCELLANEOUS INTEGER SYSTEM PROPERTIES
! LT_N      = NUMBER OF NODES PER ELEMENT
! LT_FREE   = NUMBER OF DEGREES OF FREEDOM PER ELEMENT
! LT_GEOM   = NUMBER OF GEOMETRY NODES
! N_MAT     = NUMBER OF MATERIAL TYPES
! LT_PARM   = DIMENSION OF PARAMETRIC SPACE
! 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   = UNIT FOR POST SOLUTION MATRICES STORAGE
! PRT_L_PT  = REAL PROPERTIES AT ELEMENT NODES
! PRT_MAT   = REAL ELEM PROPERTIES BASED ON MATERIAL NUMBER
! PT        = QUADRATURE COORDINATES
! S         = ELEMENT SQUARE MATRIX
! STRAIN    = STRAIN OR GRADIENT VECTOR (PLUS TWO EXTRA)
! STRAIN_0  = INITIAL STRAIN OR GRADIENT VECTOR
! STRESS    = STRESS VECTOR (PLUS TWO EXTRA)
! V         = VECTOR INTERPOLATION FUNCTIONS
! WT        = QUADRATURE WEIGHTS
! XYZ       = SPACE COORDINATES AT A POINT
! .....................................................
!   ***  ELEM_COL_MATRIX PROBLEM DEPENDENT STATEMENTS FOLLOW ***
! .....................................................

  include 'my_el_col_inc'  ! for optional C array

!     Allow a CONTAINS statement (no statements after include)
! IF ( DEBUG_EL_COL ) PRINT *, 'Leaving EL_COL ', THIS_EL
! suppress compiler warnings (never true)
! IF ( TOUCH ) CALL TOUCH_UNUSED_ITEMS_1 (AJ, AJ_INV, B, BODY, DGH, &
!              DGV, E, EB, G, H_INTG, IE, L_PT_PROP, PRT_L_PT,      &
!              PRT_MAT, STRAIN, STRAIN_0, STRESS, XYZ) 
END SUBROUTINE ELEM_COL_MATRIX
                              
SUBROUTINE SEG_COL_MATRIX (E, H_INTG, PRT_L_PT, PRT_MAT,&
                           L_PT_PROP, IE, FLUX)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * *
!          GENERATE FLUX SEGMENT COLUMN MATRIX
! * * * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants
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 Sys_Properties_Data  ! for FLT_MISC, MISC_FIX
Use Geometric_Properties 
Use Interface_Header     ! for functions
Use Select_source        ! for analytic data
 IMPLICIT NONE

 INTEGER,  INTENT(IN) :: IE

!    OPTIONAL PROPERTY VALUES
 REAL(DP), INTENT(INOUT) :: E (N_R_B, N_R_B)
 REAL(DP), INTENT(IN)    :: PRT_L_PT (LT_N, N_NP_FLO), &
                            PRT_MAT (MISC_FL) 
 INTEGER,  INTENT(IN)    :: L_PT_PROP (LT_N, N_NP_FIX)

!    OPTIONAL FLUX VALUES IF BOUNDARY SEGMENT ELEMENT
 REAL(DP), INTENT(IN) :: FLUX (L_B_N, N_G_FLUX)

!         Automatic Arrays
 REAL(DP) :: AJ (N_SPACE, N_SPACE), AJ_INV (N_SPACE, N_SPACE), &
             XYZ (N_SPACE), BODY (N_SPACE)
 REAL(DP) :: STRAIN_0 (N_R_B),          STRAIN (N_R_B + 2),      &
             STRESS (N_R_B + 2)
 REAL(DP) :: B (N_R_B, LT_FREE), EB (N_R_B, LT_FREE), H_INTG (LT_N), &
             DGH (N_SPACE, LT_N), DGV (N_SPACE, LT_FREE)
 
!                 VARIABLES:
! AJ        = JACOBIAN
! AJ_INV    = JACOBIAN INVERSE
! B         = STRAIN-DISPLACEMENT (GRADIENT) MATRIX
! BODY      = BODY FORCE VECTOR
! COORD     = SPATIAL COORDINATES OF ELEMENT'S NODES
! D         = NODAL PARAMETERS ASSOCIATED WITH AN ELEMENT
! 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
! E         = CONSTITUTIVE MATRIX
! EB        = PRODUCT OF E*B
! FLT_MISC  = SYSTEM STORAGE OF FLOATING PT MISC PROP
! G         = GEOMETRIC INTERPOLATION FUNCTIONS
! H         = SCALAR INTERPOLATION FUNCTIONS
! H_INTG    = INTEGRAL OF INTERPOLATION FUNCTIONS
! L_PT_PROP = INTEGER PROPERTIES AT EACH ELEMENT NODE
! MISC_FIX  = MISCELLANEOUS INTEGER SYSTEM PROPERTIES
! LT_N      = NUMBER OF NODES PER ELEMENT
! LT_FREE   = NUMBER OF DEGREES OF FREEDOM PER ELEMENT
! LT_GEOM   = NUMBER OF GEOMETRY NODES
! N_MAT     = NUMBER OF MATERIAL TYPES
! LT_PARM   = DIMENSION OF PARAMETRIC SPACE
! 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   = UNIT FOR POST SOLUTION MATRICES STORAGE
! PRT_L_PT  = REAL PROPERTIES AT ELEMENT NODES
! PRT_MAT   = REAL ELEM PROPERTIES BASED ON MATERIAL NUMBER
! PT        = QUADRATURE COORDINATES
! S         = ELEMENT SQUARE MATRIX
! STRAIN    = STRAIN OR GRADIENT VECTOR (PLUS TWO EXTRA)
! STRAIN_0  = INITIAL STRAIN OR GRADIENT VECTOR
! STRESS    = STRESS VECTOR (PLUS TWO EXTRA)
! V         = VECTOR INTERPOLATION FUNCTIONS
! WT        = QUADRATURE WEIGHTS
! XYZ       = SPACE COORDINATES AT A POINT
! .....................................................
!   ***  SEG_COL_MATRIX PROBLEM DEPENDENT STATEMENTS FOLLOW ***
! .....................................................

  include 'my_seg_col_inc' ! for optional C array

!     Allow a CONTAINS statement (no statements after include)
! suppress compiler warnings (never true)
! IF ( TOUCH ) CALL TOUCH_UNUSED_ITEMS_1 (AJ, AJ_INV, B, BODY, DGH, &
!              DGV, E, EB, G, H_INTG, IE, L_PT_PROP, PRT_L_PT,      &
!              PRT_MAT, STRAIN, STRAIN_0, STRESS, XYZ) 
! IF ( TOUCH ) PRINT *, FLUX
END SUBROUTINE SEG_COL_MATRIX
                              
SUBROUTINE ELEM_POST_DATA (E, H_INTG, PRT_L_PT, PRT_MAT, &
                           L_PT_PROP, IE, FLUX)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
! GENERATE DATA FOR ELEMENT POST-SOLUTION USE IN POST_PROCESS_ELEM
!           IF THAT WAS NOT DONE IN ELEM_SQ_MATRIX
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants
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 Sys_Properties_Data  ! for FLT_MISC, MISC_FIX
Use Interface_Header     ! for functions
Use Select_source        ! for analytic data
 IMPLICIT NONE

 INTEGER,  INTENT(IN)    :: IE

!    OPTIONAL PROPERTY VALUES
 REAL(DP), INTENT(INOUT) :: E (N_R_B, N_R_B)
 REAL(DP), INTENT(IN)    :: PRT_L_PT (LT_N, N_NP_FLO),  &
                            PRT_MAT (MISC_FL) 
 INTEGER,  INTENT(IN)    :: L_PT_PROP (LT_N, N_NP_FIX)

!    OPTIONAL FLUX VALUES IF BOUNDARY SEGMENT ELEMENT
 REAL(DP), INTENT(IN) :: FLUX (L_B_N, N_G_FLUX)

!         Automatic Arrays
 REAL(DP) :: AJ (N_SPACE, N_SPACE), AJ_INV (N_SPACE, N_SPACE),   &
             XYZ (N_SPACE)
 REAL(DP) :: H_INTG (LT_N   ),           &
             B (N_R_B, LT_FREE),     EB (N_R_B, LT_FREE),        &
             STRAIN_0 (N_R_B),       STRAIN (N_R_B + 2),         &
             STRESS (N_R_B + 2),     BODY (N_SPACE)
 REAL(DP) :: DGH (N_SPACE, LT_N), DGV (N_SPACE, LT_FREE)

!                 VARIABLES:
! B        = STRAIN-DISPLACEMENT (GRADIENT) MATRIX
! BODY     = BODY FORCE VECTOR
! D        = NODAL PARAMETERS ASSOCIATED WITH AN ELEMENT
! DGH      = GLOBAL DERIVATIVES SCALAR INTERPOLATION FUNCTIONS
! DGV       = GLOBAL DERIVATIVES VECTOR INTERPOLATION FUNCTIONS
! E        = CONSTITUTIVE MATRIX
! EB       = PRODUCT OF E*B
! FLT_MISC = SYSTEM STORAGE OF FLOATING PT MISC PROP
! H        = SCALAR INTERPOLATION FUNCTIONS
! H_INTG   = INTEGRAL OF INTERPOLATION FUNCTIONS
! L_PT_PROP= INTEGER PROPERTIES AT EACH ELEMENT NODE
! MISC_FIX = MISCELLANEOUS INTEGER SYSTEM PROPERTIES
! LT_N     = NUMBER OF NODES PER ELEMENT
! LT_FREE  = NUMBER OF DEGREES OF FREEDOM PER ELEMENT
! LT_GEOM  = NUMBER OF GEOMETRY NODES
! N_MAT    = NUMBER OF MATERIAL TYPES
! LT_PARM  = DIMENSION OF PARAMETRIC SPACE
! 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  = UNIT FOR POST SOLUTION MATRICES STORAGE
! PRT_L_PT = REAL PROPERTIES AT ELEMENT NODES
! PRT_MAT  = REAL ELEM PROPERTIES BASED ON MATERIAL NUMBER
! STRAIN   = STRAIN OR GRADIENT VECTOR (PLUS TWO EXTRA)
! STRAIN_0 = INITIAL STRAIN OR GRADIENT VECTOR
! STRESS   = STRESS VECTOR (PLUS TWO EXTRA)
! V         = VECTOR INTERPOLATION FUNCTIONS
! .....................................................
!    *** ELEM_POST_DATA PROBLEM DEPENDENT STATEMENTS FOLLOW ***
! .....................................................

  include 'my_el_post_inc' ! is rarely needed

! suppress compiler warnings (never true)
  IF ( TOUCH ) CALL TOUCH_UNUSED_ITEMS_1 (AJ, AJ_INV, B, BODY, DGH, &
               DGV, E, EB, G, H_INTG, IE, L_PT_PROP, PRT_L_PT,      &
               PRT_MAT, STRAIN, STRAIN_0, STRESS, XYZ) 
  IF ( TOUCH ) PRINT *, FLUX
END SUBROUTINE ELEM_POST_DATA
                             
SUBROUTINE POST_PROCESS_ELEM (PRT_L_PT, PRT_MAT, &
                              L_PT_PROP, ITER, IE,&
                              DD, DD_OLD) !b
! * * * * * * * * * * * * * * * * * * * * * * * * * * * *
! ELEMENT LEVEL POST-SOLUTION CALCULATIONS, FROM ELEM_POST_DATA
! * * * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants
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 Sys_Properties_Data  ! for FLT_MISC, MISC_FIX
Use Interface_Header     ! for functions
Use Select_source        ! for analytic data
 IMPLICIT NONE
 INTEGER,  INTENT(IN) :: ITER, IE ! ITERATION, ELEMENT NUMBER
 REAL(DP), INTENT(INOUT) :: DD (N_D_FRE), DD_OLD (N_D_FRE)

!    OPTIONAL PROPERTY AND SOLUTION VALUES
 REAL(DP), INTENT(IN) :: PRT_L_PT (LT_N, N_NP_FLO), &
                         PRT_MAT (MISC_FL) 
 INTEGER,  INTENT(IN) :: L_PT_PROP (LT_N, N_NP_FIX)

!         Automatic Arrays (Usually used.)
 REAL(DP) :: B (N_R_B, LT_FREE), DGH (N_SPACE, LT_N),          &
             DGV (N_SPACE, LT_FREE),                           &
             E (N_R_B, N_R_B), EB (N_R_B, LT_FREE),            &
             STRAIN (N_R_B + 2), STRAIN_0 (N_R_B),             &
             STRESS (N_R_B + 2), BODY (N_SPACE)

!    OPTIONAL FOR NUMERICAL INTEGRATION
 REAL(DP) :: AJ (N_SPACE, N_SPACE), AJ_INV (N_SPACE, N_SPACE), &
             H_INTG (LT_N),         XYZ (N_SPACE)

!                 VARIABLES:
! AJ        = JACOBIAN
! AJ_INV    = JACOBIAN INVERSE
! B         = STRAIN-DISPLACEMENT (GRADIENT) MATRIX
! BODY      = BODY FORCE VECTOR
! COORD     = SPATIAL COORDINATES OF ELEMENT'S NODES
! D         = NODAL PARAMETERS ASSOCIATED WITH AN ELEMENT
! 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
! E         = CONSTITUTIVE MATRIX
! EB        = PRODUCT OF E*B
! FLT_MISC  = SYSTEM STORAGE OF FLOATING PT MISC PROP
! G         = GEOMETRIC INTERPOLATION FUNCTIONS
! H         = SCALAR INTERPOLATION FUNCTIONS
! H_INTG    = INTEGRAL OF INTERPOLATION FUNCTIONS
! ITER      = CURRENT ITERATION NUMBER
! L_PT_PROP = INTEGER PROPERTIES AT EACH ELEMENT NODE
! MISC_FIX  = MISCELLANEOUS INTEGER SYSTEM PROPERTIES
! LT_N      = NUMBER OF NODES PER ELEMENT
! LT_FREE   = NUMBER OF DEGREES OF FREEDOM PER ELEMENT
! LT_GEOM   = NUMBER OF GEOMETRY NODES
! N_ITER    = MAX NUMBER OF ITERATIONS
! N_MAT     = NUMBER OF MATERIAL TYPES
! N_MIXED   = NUMBER OF SEGMENTS WITH MIXED BOUNDARY CONDITIONS
! LT_PARM   = DIMENSION OF PARAMETRIC SPACE
! 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   = UNIT FOR POST SOLUTION MATRICES STORAGE
! PRT_L_PT  = REAL PROPERTIES AT ELEMENT NODES
! PRT_MAT   = REAL ELEM PROPERTIES BASED ON MATERIAL NUMBER
! PT        = QUADRATURE COORDINATES
! S         = ELEMENT SQUARE MATRIX
! STRAIN    = STRAIN OR GRADIENT VECTOR (PLUS TWO EXTRA)
! STRAIN_0  = INITIAL STRAIN OR GRADIENT VECTOR
! STRESS    = STRESS VECTOR (PLUS TWO EXTRA)
! V         = VECTOR INTERPOLATION FUNCTIONS
! WT        = QUADRATURE WEIGHTS
! XYZ       = SPACE COORDINATES AT A POINT
! ....................................................
!  *** POST_PROCESS_ELEM PROBLEM DEPENDENT STATEMENTS FOLLOW ***
! ....................................................

  include 'my_post_el_inc'  ! is usually needed

!! suppress compiler warnings (touch is never true)
!  IF ( TOUCH ) CALL TOUCH_UNUSED_ITEMS_1 (AJ, AJ_INV, B, BODY, DGH, &
!               DGV, E, EB, G, H_INTG, IE, L_PT_PROP, PRT_L_PT,      &
!               PRT_MAT, STRAIN, STRAIN_0, STRESS, XYZ)
!  IF ( TOUCH ) PRINT *, ITER
END SUBROUTINE POST_PROCESS_ELEM
                                
SUBROUTINE POST_PROCESS_MIXED (PRT_L_PT, PRT_MAT, &
                               L_PT_PROP, ITER, IE)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * *
! MIXED SEGMENT POST-SOLUTION CALCULATIONS
! * * * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants
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 Sys_Properties_Data  ! for FLT_MISC, MISC_FIX
Use Interface_Header     ! for functions
Use Select_source        ! for analytic data
 IMPLICIT NONE
 INTEGER,  INTENT(IN) :: ITER, IE ! ITERATION, ELEMENT NUMBER

!    OPTIONAL PROPERTY AND SOLUTION VALUES
 REAL(DP), INTENT(IN) :: PRT_L_PT (LT_N, N_NP_FLO), &
                         PRT_MAT (MISC_FL) 
 INTEGER,  INTENT(IN) :: L_PT_PROP (LT_N, N_NP_FIX)

!         Automatic Arrays (Usually used.)
 REAL(DP) :: B (N_R_B, LT_FREE), DGH (N_SPACE, LT_N),          &
             DGV (N_SPACE, LT_FREE),                           &
             E (N_R_B, N_R_B), EB (N_R_B, LT_FREE),            &
             STRAIN (N_R_B + 2), STRAIN_0 (N_R_B),             &
             STRESS (N_R_B + 2), BODY (N_SPACE)

!    OPTIONAL FOR NUMERICAL INTEGRATION
 REAL(DP) :: AJ (N_SPACE, N_SPACE), AJ_INV (N_SPACE, N_SPACE), &
             H_INTG (LT_N),         XYZ (N_SPACE)

!                 VARIABLES:
! AJ        = JACOBIAN
! AJ_INV    = JACOBIAN INVERSE
! B         = STRAIN-DISPLACEMENT (GRADIENT) MATRIX
! BODY      = BODY FORCE VECTOR
! COORD     = SPATIAL COORDINATES OF ELEMENT'S NODES
! D         = NODAL PARAMETERS ASSOCIATED WITH AN ELEMENT
! 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
! E         = CONSTITUTIVE MATRIX
! EB        = PRODUCT OF E*B
! FLT_MISC  = SYSTEM STORAGE OF FLOATING PT MISC PROP
! G         = GEOMETRIC INTERPOLATION FUNCTIONS
! H         = SCALAR INTERPOLATION FUNCTIONS
! H_INTG    = INTEGRAL OF INTERPOLATION FUNCTIONS
! ITER      = CURRENT ITERATION NUMBER
! L_PT_PROP = INTEGER PROPERTIES AT EACH ELEMENT NODE
! MISC_FIX  = MISCELLANEOUS INTEGER SYSTEM PROPERTIES
! LT_N      = NUMBER OF NODES PER ELEMENT
! LT_FREE   = NUMBER OF DEGREES OF FREEDOM PER ELEMENT
! LT_GEOM   = NUMBER OF GEOMETRY NODES
! N_ITER    = MAX NUMBER OF ITERATIONS
! N_MAT     = NUMBER OF MATERIAL TYPES
! N_MIXED   = NUMBER OF SEGMENTS WITH MIXED BOUNDARY CONDITIONS
! LT_PARM   = DIMENSION OF PARAMETRIC SPACE
! 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   = UNIT FOR POST SOLUTION MATRICES STORAGE
! PRT_L_PT  = REAL PROPERTIES AT ELEMENT NODES
! PRT_MAT   = REAL ELEM PROPERTIES BASED ON MATERIAL NUMBER
! PT        = QUADRATURE COORDINATES
! S         = ELEMENT SQUARE MATRIX
! STRAIN    = STRAIN OR GRADIENT VECTOR (PLUS TWO EXTRA)
! STRAIN_0  = INITIAL STRAIN OR GRADIENT VECTOR
! STRESS    = STRESS VECTOR (PLUS TWO EXTRA)
! V         = VECTOR INTERPOLATION FUNCTIONS
! WT        = QUADRATURE WEIGHTS
! XYZ       = SPACE COORDINATES AT A POINT
! ....................................................
!  *** POST_PROCESS_MIXED PROBLEM DEPENDENT STATEMENTS FOLLOW ***
! ....................................................

  include 'my_post_mixed_inc'  ! is usually needed

!! suppress compiler warnings (touch is never true)
!  IF ( TOUCH ) CALL TOUCH_UNUSED_ITEMS_1 (AJ, AJ_INV, B, BODY, DGH, &
!               DGV, E, EB, G, H_INTG, IE, L_PT_PROP, PRT_L_PT,      &
!               PRT_MAT, STRAIN, STRAIN_0, STRESS, XYZ)
!  IF ( TOUCH ) PRINT *, ITER
END SUBROUTINE POST_PROCESS_MIXED
                                
FUNCTION START_DOF_VALUE (IG, XYZ)  
! * * * * * * * * * * * * * * * * * * * * * * * * * *
! DEFINE STARTING VALUE OF PARAMETER IG IN TERMS OF
! COORDINATES OF THE NODE  (FOR ITERATIVE SOLUTIONS)
! * * * * * * * * * * * * * * * * * * * * * * * * * *
!    A PROBLEM DEPENDENT ROUTINE
Use System_Constants
Use Select_source       ! for analytic data
Use Sys_Properties_Data ! for GET_INTEGER_* or GET_REAL_* functions
IMPLICIT NONE
 INTEGER,  INTENT(IN) :: IG               ! local dof number
 REAL(DP), INTENT(IN) :: XYZ (N_SPACE)    ! position
 REAL(DP)             :: START_DOF_VALUE  ! result        

! IG      = LOCAL PARAMETER NUMBER AT NODE
! N_SPACE = DIMENSION OF SPACE
! XYZ     = SPATIAL COORDINATE ARRAY OF NODE
! ....................................................
! ** PROBLEM DEPENDENT START STATEMENTS FOLLOW **
! ....................................................

  include 'my_iter_start_inc'

! suppress compiler warnings (touch is never true)
  IF ( TOUCH ) PRINT *,'START_DOF_VALUE:', IG, XYZ
END FUNCTION START_DOF_VALUE 
                  
FUNCTION GET_APPLICATION_B_MATRIX (DGH, XYZ) RESULT (B)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
!  GET APPLICATION DEPENDENT B MATRIX FOR STRAIN OR GRADIENTS
!    (USED IF SUPERCONVERGENT PATCH GRADIENTS ARE ACTIVE)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants
Use Elem_Type_Data  ! for LT_N, LT_FREE, & H (LT_N) IF NEEDED
! Use Interface_Header     ! for functions
 IMPLICIT NONE
 REAL (DP), INTENT(IN)  :: DGH (N_SPACE, LT_N)
 REAL (DP), INTENT(IN)  :: XYZ (N_SPACE)
 REAL (DP)              :: B   (N_R_B,   LT_FREE)

!   B          = GRADIENT VERSUS DOF MATRIX (N_R_B, LT_FREE)
!   DGH        = GLOBAL DERIVS OF SCALAR FUNCTIONS H (N_SPACE, LT_N)
!   H          = INTERPOLATION FUNCTIONS FOR AN ELEMENT SCALAR
!_t LT_FREE    = NUMBER OF DEGREES OF FREEDOM PER ELEMENT
!_t LT_N       = MAXIMUM NUMBER OF NODES FOR ELEMENT TYPE
!_s N_R_B      = NUMBER OF ROWS IN B AND E MATRICES
!_s N_SPACE    = DIMENSION OF SPACE
!   XYZ        = SPACE COORDINATES AT A POINT (N_SPACE)

! INTERFACE ! for Poisson problems
!    FUNCTION COPY_DGH_INTO_B_MATRIX (DGH) RESULT (B)
!    Use System_Constants
!    Use Elem_Type_Data  ! for LT_N, LT_FREE  
!     REAL (DP), INTENT(IN)  :: DGH (N_SPACE, LT_N)
!     REAL (DP)              :: B   (N_R_B,   LT_FREE)
!    END FUNCTION COPY_DGH_INTO_B_MATRIX
! END INTERFACE

   include 'my_b_matrix_inc'  !if application specific, or
!                  or
!  CALL B_MATRIX_EQUAL_DGH (DGH, B) ! no INTERFACE required
! suppress compiler warnings (touch is never true)
  IF ( TOUCH ) PRINT *,DGH (1, 1), XYZ (1), B (1, 1)
  !b STOP 'Edit GET_APPLICATION_B_MATRIX to use my_b_matrix_inc'
END FUNCTION GET_APPLICATION_B_MATRIX 

SUBROUTINE APPLICATION_B_MATRIX (DGH, XYZ, B)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
!  GET APPLICATION DEPENDENT B MATRIX FOR STRAIN OR GRADIENTS
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants
Use Elem_Type_Data  ! for LT_N, LT_FREE, H (LT_N), and     
                    ! COORD (LT_N, N_SPACE) if needed
Use Interface_Header     ! for functions
 IMPLICIT NONE
 REAL (DP), INTENT(IN)  :: DGH (N_SPACE, LT_N)
 REAL (DP), INTENT(IN)  :: XYZ (N_SPACE)
 REAL (DP), INTENT(OUT) :: B   (N_R_B,   LT_FREE)

!b INTEGER :: I, J, K, L, M, N  ! for loops

!   B       = GRADIENT VERSUS DOF MATRIX (N_R_B, LT_FREE)
!   COORD   = COORDINATES OF ALL NODES ON ELEMENT (LT_N, N_SPACE)
!   DGH     = GLOBAL DERIVS OF SCALAR FUNCTIONS H (N_SPACE, LT_N)
!_t LT_FREE = NUMBER OF DEGREES OF FREEDOM PER ELEMENT
!_t LT_N    = MAXIMUM NUMBER OF NODES FOR ELEMENT TYPE
!_s N_R_B   = NUMBER OF ROWS IN B AND E MATRICES
!_s N_SPACE = DIMENSION OF SPACE
!   XYZ     = SPACE COORDINATES AT A POINT (N_SPACE)

! INTERFACE ! for Poisson problems
!    FUNCTION COPY_DGH_INTO_B_MATRIX (DGH) RESULT (B)
!    Use System_Constants
!    Use Elem_Type_Data  ! for LT_N, LT_FREE  
!     REAL (DP), INTENT(IN)  :: DGH (N_SPACE, LT_N)
!     REAL (DP)              :: B   (N_R_B,   LT_FREE)
!    END FUNCTION COPY_DGH_INTO_B_MATRIX
! END INTERFACE

   include 'my_b_matrix_inc'  ! application specific

! suppress compiler warnings (touch is never true)
  IF ( TOUCH ) PRINT *, DGH (1, 1), XYZ (1), B (1, 1)
  !b STOP 'Edit APPLICATION_B_MATRIX to use my_b_matrix_inc'
END SUBROUTINE APPLICATION_B_MATRIX 

FUNCTION GET_APPLICATION_E_MATRIX (IE, XYZ) RESULT (E)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
!  GET APPLICATION DEPENDENT E MATRIX FOR CONSTITUTIVE LAW   
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants
Use Sys_Properties_Data
Use Elem_Type_Data  ! for D (LT_FREE), the nodal solution values
 IMPLICIT NONE
 INTEGER,  INTENT (IN)  :: IE  ! or THIS_EL from System_Constants
 REAL (DP), INTENT(IN)  :: XYZ (N_SPACE)
 REAL (DP)              :: E   (N_R_B, N_R_B)     

!b INTEGER :: I, J, K, L, M, N  ! for loops

!   D          = NODAL PARAMETERS ASSOCIATED WITH AN ELEMENT
!   E          = CONSTITUTIVE MATRIX
!   IE         = ELEMENT NUMBER
!_s N_R_B      = NUMBER OF ROWS IN B AND E MATRICES
!_s N_SPACE    = DIMENSION OF SPACE
!   XYZ        = SPACE COORDINATES AT A POINT (N_SPACE)

! INTERFACE ! for Poisson or constant arrays
!    FUNCTION GET_REAL_IDENTITY (N)  RESULT (EYE)
!    Use System_Constants
!     INTEGER,  INTENT(IN)  :: N            ! SIZE OF MATRIX
!     REAL(DP)              :: EYE (N, N)   ! IDENTITY MATRIX
!    END FUNCTION GET_REAL_IDENTITY
! END INTERFACE

   include 'my_e_matrix_inc'

!  CALL REAL_IDENTITY (N_R_B, E) ! DEFAULT TO IDENTITY MATRIX or
!  STOP 'Edit GET_APPLICATION_E_MATRIX to use my_e_matrix_inc'

! suppress compiler warnings (touch is never true)
  IF ( TOUCH ) PRINT *, IE, XYZ (1), E (1, 1)
END FUNCTION GET_APPLICATION_E_MATRIX 

SUBROUTINE APPLICATION_E_MATRIX (IE, XYZ, E)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
!  GET APPLICATION DEPENDENT E MATRIX FOR CONSTITUTIVE LAW   
!    (USED IF SUPERCONVERGENT PATCH GRADIENTS ARE ACTIVE)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants
Use Sys_Properties_Data
Use Elem_Type_Data  ! for D (LT_FREE), the nodal solution values
 IMPLICIT NONE
 INTEGER,  INTENT (IN)  :: IE
 REAL (DP), INTENT(IN)  :: XYZ (N_SPACE)
 REAL (DP), INTENT(OUT) :: E   (N_R_B, N_R_B)     

!b INTEGER :: I, J, K, L, M, N  ! for loops

!   D          = NODAL PARAMETERS ASSOCIATED WITH AN ELEMENT
!   E          = CONSTITUTIVE MATRIX
!   IE         = ELEMENT NUMBER
!_s N_R_B      = NUMBER OF ROWS IN B AND E MATRICES
!_s N_SPACE    = DIMENSION OF SPACE
!   XYZ        = SPACE COORDINATES AT A POINT (N_SPACE)

   include 'my_e_matrix_inc'
!  CALL REAL_IDENTITY (N_R_B, E) ! DEFAULT TO IDENTITY MATRIX or
!  STOP 'Edit APPLICATION_E_MATRIX to use my_e_matrix_inc'

! suppress compiler warnings (touch is never true)
  IF ( TOUCH ) PRINT *, IE, XYZ (1), E (1, 1)
END SUBROUTINE APPLICATION_E_MATRIX

SUBROUTINE E_MATRIX_FROM_XYZ (XYZ, E)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
!  GET APPLICATION DEPENDENT E MATRIX FOR CONSTITUTIVE LAW   
!          FROM THE SPATIAL COORDINATES AT A POINT
!    (USED IF SUPERCONVERGENT PATCH GRADIENTS ARE ACTIVE)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants
Use Sys_Properties_Data
Use Elem_Type_Data  ! for D (LT_FREE), the nodal solution values
 IMPLICIT NONE
 REAL (DP), INTENT(IN)  :: XYZ (N_SPACE)
 REAL (DP), INTENT(OUT) :: E   (N_R_B, N_R_B)     

!   D          = NODAL PARAMETERS ASSOCIATED WITH AN ELEMENT
!   E          = CONSTITUTIVE MATRIX
!_s N_R_B      = NUMBER OF ROWS IN B AND E MATRICES
!_s N_SPACE    = DIMENSION OF SPACE
!   XYZ        = SPACE COORDINATES AT A POINT (N_SPACE)
   include 'my_e_matrix_inc'
   !b IF ( DEBUG_E ) PRINT *,'Entered E_MATRIX_FROM_DOF'

!  STOP 'Edit E_MATRIX_FROM_XYZ to use my_e_matrix_inc'

! suppress compiler warnings (touch is never true)
  IF ( TOUCH ) PRINT *, XYZ (1), E (1, 1)
END SUBROUTINE E_MATRIX_FROM_XYZ

SUBROUTINE E_MATRIX_FROM_DOF (DOF, E)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
!  GET APPLICATION DEPENDENT E MATRIX FOR CONSTITUTIVE LAW   
!          FROM THE DOF INTERPOLATED AT A POINT
!    (USED IF SUPERCONVERGENT PATCH GRADIENTS ARE ACTIVE)
!           See INTERPOLATE_DOF_WITH_H
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants
Use Sys_Properties_Data
Use Elem_Type_Data  ! for D (LT_FREE), the nodal solution values
 IMPLICIT NONE
 REAL (DP), INTENT(IN)  :: DOF (N_G_DOF)
 REAL (DP), INTENT(OUT) :: E   (N_R_B, N_R_B)     

!   D          = NODAL PARAMETERS ASSOCIATED WITH AN ELEMENT
!   E          = CONSTITUTIVE MATRIX
!_t LT_FREE    = DEGREES OF FREEDOM PER ELEMENT, N_G_DOF*LT_N
!_t LT_N       = MAXIMUM NUMBER OF NODES FOR ELEMENT TYPE
!_s N_G_DOF    = NUMBER OF GENERALIZED PARAMETERS (DOF) PER NODE
!_s N_R_B      = NUMBER OF ROWS IN B AND E MATRICES
!   DOF        = INTERPOLATED VALUE, H * D              
!   include 'my_e_matrix_inc'
   !b IF ( DEBUG_E ) PRINT *,'Entered E_MATRIX_FROM_DOF'

  STOP 'Edit E_MATRIX_FROM_DOF to use my_e_matrix_inc'

! suppress compiler warnings (touch is never true)
  IF ( TOUCH ) PRINT *, DOF (1), E (1, 1)
END SUBROUTINE E_MATRIX_FROM_DOF

!SUBROUTINE GET_CONSTANT_E (IE, E) !b turned off 4/8/02
!! * * * * * * * * * * * * * * * * * * * * * * * * * * * *
!!      COMPUTE A CONSTANT CONSTITUTIVE MATRIX, E
!! * * * * * * * * * * * * * * * * * * * * * * * * * * * *
!Use System_Constants
!Use Elem_Type_Data  ! for: LT_N, LT_PARM, COORD (LT_N, N_SPACE), &
!                    ! DLH (LT_PARM, LT_N),  H (LT_N), etc
!Use Sys_Properties_Data  ! for FLT_MISC, MISC_FIX
!!b Use Interface_Header
! IMPLICIT NONE
!
! INTEGER,  INTENT(IN)    :: IE                 ! ELEMENT NUMBER
! REAL(DP), INTENT(INOUT) :: E (N_R_B, N_R_B)   ! CONSTITUTIVE MATRIX
!
!! INTERFACE ! for Poisson or constant arrays
!!    FUNCTION GET_REAL_IDENTITY (N)  RESULT (EYE)
!!    Use System_Constants
!!     INTEGER,  INTENT(IN)  :: N            ! SIZE OF MATRIX
!!     REAL(DP)              :: EYE (N, N)   ! IDENTITY MATRIX
!!    END FUNCTION GET_REAL_IDENTITY
!! END INTERFACE
!
!    include 'my_e_matrix_inc'
!
!! suppress compiler warnings (touch is never true)
!  IF ( TOUCH ) PRINT *, IE, E (1, 1)
!  STOP 'Edit GET_CONSTANT_E to use my_e_matrix_inc'
!END SUBROUTINE GET_CONSTANT_E

SUBROUTINE GET_E_AT_QP (IE, IQ, VALUES, E)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * *
! COMPUTE CONSTITUTIVE MATRIX, E, AT A QUADRATURE POINT
! * * * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants
Use Elem_Type_Data  ! for: LT_N, LT_PARM, COORD (LT_N, N_SPACE), &
                    ! DLH (LT_PARM, LT_N),  H (LT_N), D (LT_FREE)
Use Sys_Properties_Data  ! for FLT_MISC, MISC_FIX
Use Interface_Header
 IMPLICIT NONE

 INTEGER,  INTENT(IN)    :: IE, IQ                 ! ELEM, POINT
 REAL(DP), INTENT(IN)    :: VALUES (N_NP_FLO)      ! PROPERTIES AT PT
 REAL(DP), INTENT(INOUT) :: E      (N_R_B, N_R_B)  ! CONSTITUTIVE

! OVERWRITE WITH APPLICATION FORM USING VALUES
! include 'my_e_matrix_inc'
! suppress compiler warnings (touch is never true)
  IF ( TOUCH ) PRINT *, IE, IQ, VALUES (1), E (1, 1)
  STOP 'Edit GET_E_AT_QP to use my_e_matrix_inc'
END SUBROUTINE GET_E_AT_QP

SUBROUTINE GET_E_AT_XYZ (XYZ, E)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * *
! COMPUTE CONSTITUTIVE MATRIX, E, AT A PHYSICAL POINT
! * * * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants
Use Elem_Type_Data  ! for: LT_N, LT_PARM, COORD (LT_N, N_SPACE), &
                    ! DLH (LT_PARM, LT_N),  H (LT_N), D (LT_FREE)
Use Sys_Properties_Data  ! for FLT_MISC, MISC_FIX
Use Interface_Header
 IMPLICIT NONE

 REAL(DP), INTENT(IN)    :: XYZ (N_SPACE)       ! LOCAL COORDINATES  
 REAL(DP), INTENT(INOUT) :: E   (N_R_B, N_R_B)  ! CONSTITUTIVE

! OVERWRITE WITH APPLICATION FORM USING VALUES
! include 'my_e_matrix_at_xyz_inc'
! suppress compiler warnings (touch is never true)
  IF ( TOUCH ) PRINT *, XYZ
  STOP 'Edit GET_E_AT_XYZ to use my_e_matrix_at_xyz_inc'   
END SUBROUTINE GET_E_AT_XYZ

SUBROUTINE TOUCH_UNUSED_ITEMS_1 (AJ, AJ_INV, B, BODY, DGH, DGV, &
           E, EB, G_FAKE, H_INTG, IE, L_PT_PROP, PRT_L_PT, PRT_MAT,  &
           STRAIN, STRAIN_0, STRESS, XYZ)
! ..............................................................
! fake touch of unused items to stop compiler warnings
! ..............................................................
Use System_Constants
Use Elem_Type_Data  ! for: LT_N, LT_GEOM, G(LT_GEOM)
 IMPLICIT NONE
 REAL(DP), INTENT(IN) :: AJ (N_SPACE, N_SPACE),                      &
   AJ_INV (N_SPACE, N_SPACE), B (N_R_B, N_EL_FRE), BODY (N_SPACE),   &
   DGH (N_SPACE, NOD_PER_EL), DGV (N_SPACE, N_EL_FRE),               &
   E (N_R_B, N_R_B), EB (N_R_B, N_EL_FRE), G_FAKE (N_GEOM),          &
   H_INTG (NOD_PER_EL), PRT_L_PT (NOD_PER_EL, N_NP_FLO),             &
   PRT_MAT (MISC_FL), STRAIN (N_R_B + 2),  STRAIN_0 (N_R_B),         &
   STRESS (N_R_B + 2), XYZ (N_SPACE)
 INTEGER,  INTENT(IN) :: IE, L_PT_PROP (LT_N, N_NP_FIX)

   IF ( TOUCH ) THEN ! THIS IS NEVER TRUE
     PRINT *,'WARNING: TOUCH_UNUSED_ITEMS_1, TOUCH IS NEVER TRUE'
     N_WARN = N_WARN + 1
     PRINT *, AJ (1, 1), AJ_INV (1, 1), B (1, 1), BODY (1),          &
       DGH (1, 1), DGV (1, 1), E (1, 1), EB (1, 1), G_FAKE (1),      &
       H_INTG (1), PRT_L_PT (1, 1), PRT_MAT (1), STRAIN (1),         &
       STRAIN_0 (1), STRESS (1), XYZ (1), IE, L_PT_PROP (1,1)
   END IF ! NEVER DONE
           
END SUBROUTINE TOUCH_UNUSED_ITEMS_1

! End application_lib.f