! copyright 2005,  J. E. Akin, all rights reserved.

SUBROUTINE B_MATRIX_AXISYM_ELASTIC (NOD_PER_EL, N_SPACE,     &
                                    N_G_DOF, GDH, H, R, NS, B)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
! AXISYMMETRIC ELASTICITY STRAIN-DISPLACEMENT RELATIONS  (B)
! STRESS & STRAIN COMPONENT ORDER: RR, ZZ, RZ, AND TT      
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
Use Precision_Module
 IMPLICIT NONE
 INTEGER,  INTENT(IN)  :: NOD_PER_EL, N_SPACE, N_G_DOF, NS
 REAL(DP), INTENT(IN)  :: R, GDH (N_SPACE, NOD_PER_EL * N_G_DOF)
 REAL(DP), INTENT(IN)  :: H (NOD_PER_EL * N_G_DOF)
 REAL(DP), INTENT(OUT) :: B (NS, NOD_PER_EL * N_G_DOF)
 INTEGER :: J, K, L
                               
! B          = STRAIN-DISPLACEMENT MATRIX (RETURNED)
! GDH        = GLOBAL DERIVATIVES OF H
! H          = ELEMENT INTERPOLATION FUNCTIONS
! NOD_PER_EL = NUMBER OF NODES PER ELEMENT
! N_G_DOF    = NUMBER OF PARAMETERS PER NODE
! NS         = NUMBER OF STRAINS (ROWS IN B): XX, YY, XY, HOOP
! N_SPACE    = DIMENSION OF SPACE
! R,Z,T        DENOTE RADIAL, AXIAL, CIRCUMFERENCE
 
  B = 0.0d0  
  DO J = 1, NOD_PER_EL             ! ROW NUMBER
    K = N_G_DOF * (J - 1) + 1      ! FIRST COLUMN, U
    L = K + 1                      ! SECOND COLUMN, V
                                             
    B (1, K) = GDH (1, J)          ! DU/DX  FOR XX NORMAL
    B (3, K) = GDH (2, J)          ! DU/DY  FOR XY SHEAR
    IF ( R <= 0.0d0 ) STOP 'R=0, B_MATRIX_AXISYM_ELASTIC'
    B (4, K) = H (J) / R           ! U/R HOOP, ZZ NORMAL
    B (2, L) = GDH (2, J)          ! DV/DY  FOR YY NORMAL
    B (3, L) = GDH (1, J)          ! DV/DX  FOR XY SHEAR
  END DO
END SUBROUTINE B_MATRIX_AXISYM_ELASTIC
                        
SUBROUTINE B_MATRIX_BAR_ELASTIC (NOD_PER_EL, N_SPACE, N_G_DOF, &
                                 GDH, NS, B)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
! ELASTIC BAR STRAIN-DISPLACEMENT RELATION  (B)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
Use Precision_Module
 IMPLICIT NONE
 INTEGER,  INTENT(IN)  :: NOD_PER_EL, N_SPACE, N_G_DOF, NS
 REAL(DP), INTENT(IN)  :: GDH (N_SPACE, NOD_PER_EL * N_G_DOF)
 REAL(DP), INTENT(OUT) :: B (NS, NOD_PER_EL * N_G_DOF)
 INTEGER :: J, K
                            
! B          = STRAIN-DISPLACEMENT MATRIX (RETURNED)
! GDH        = GLOBAL DERIVATIVES OF ELEM INTERPOLATION FUNCTIONS
! NOD_PER_EL = NUMBER OF NODES PER ELEMENT
! N_G_DOF    = NUMBER OF PARAMETERS PER NODE
! NS         = NUMBER OF STRAINS (ROWS IN B) = 1
! N_SPACE    = DIMENSION OF SPACE
 
  DO J = 1, NOD_PER_EL             ! ROW NUMBER
    K = N_G_DOF * (J - 1) + 1      ! FIRST COLUMN, U
    B (1, K) = GDH (1, J)          ! DU/DX  FOR XX NORMAL
  END DO
END SUBROUTINE B_MATRIX_BAR_ELASTIC


SUBROUTINE B_MATRIX_CYL_ELASTIC (NOD_PER_EL, GDH, H, R, B)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
! ELASTIC  CYLINDER STRAIN-DISPLACEMENT RELATIONS  (B)
!      STRESS & STRAIN COMPONENT ORDER: RR, AND TT      
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
Use Precision_Module
 IMPLICIT NONE
 INTEGER,  INTENT(IN)  :: NOD_PER_EL
 REAL(DP), INTENT(IN)  :: GDH (1, NOD_PER_EL)
 REAL(DP), INTENT(IN)  :: H   (NOD_PER_EL), R
 REAL(DP), INTENT(OUT) :: B   (2, NOD_PER_EL)
 INTEGER :: J
                               
! B          = STRAIN-DISPLACEMENT MATRIX (RETURNED)
! GDH        = GLOBAL DERIVATIVES OF H
! H          = ELEMENT INTERPOLATION FUNCTIONS
! NOD_PER_EL = NUMBER OF NODES PER ELEMENT
! R,T          DENOTE RADIAL, CIRCUMFERENCE
 
  B = 0.d0  
  DO J = 1, NOD_PER_EL             ! ROW NUMBER
    B (1, J) = GDH (1, J)          ! DU/DR  FOR RR NORMAL
    IF ( R <= 0.d0 ) STOP 'R=0, B_MATRIX_CYL_ELASTIC'
    B (2, J) = H (J) / R           ! U/R HOOP
  END DO
END SUBROUTINE B_MATRIX_CYL_ELASTIC
                         
SUBROUTINE B_MATRIX_ELASTIC (I_OPT, NOD_PER_EL, N_SPACE,   &
                             N_G_DOF, GDH, H, R, NS, B)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
! ELASTICITY STRAIN-DISPLACEMENT RELATIONS  (B)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
Use Precision_Module
 IMPLICIT NONE
 INTEGER,  INTENT(IN)  :: I_OPT, NOD_PER_EL, N_SPACE, N_G_DOF, NS
 REAL(DP), INTENT(IN)  :: R, GDH (N_SPACE, NOD_PER_EL * N_G_DOF)
 REAL(DP), INTENT(IN)  :: H (NOD_PER_EL * N_G_DOF)
 REAL(DP), INTENT(OUT) :: B (NS, NOD_PER_EL * N_G_DOF)
 INTEGER :: J, K, L, M
                                  
! B          = STRAIN-DISPLACEMENT MATRIX (RETURNED)
! GDH        = GLOBAL DERIVATIVES OF H
! H          = ELEMENT INTERPOLATION FUNCTIONS
! I_OPT       = ELASTICITY CLASS
!              = 1, AXIAL BAR,    N_G_DOF = 1, NS = 1
!              = 2, PLANE STRESS, N_G_DOF = 2, NS = 3
!              = 3, PLANE STRAIN, N_G_DOF = 2, NS = 3
!              = 4, AXISYMMETRIC, N_G_DOF = 2, R = RADIUS, NS = 4
!              = 5, 3-D SOLID,    N_G_DOF = 3, NS = 6
! NOD_PER_EL = NUMBER OF NODES PER ELEMENT
! N_G_DOF    = NUMBER OF PARAMETERS PER NODE
! NS         = NUMBER OF STRAINS (ROWS IN B)
!              XX, YY, XY, ZZ, XZ, YZ, 1 OR 3 OR 4 OR 6
! N_SPACE    = DIMENSION OF SPACE
 
  B = 0.d0
  IF ( I_OPT < 1 .OR. I_OPT > 5 ) STOP 'I_OPT ERROR, B_MATRIX_ELASTIC'
  DO J = 1, NOD_PER_EL                        ! ROW NUBMER
    K = N_G_DOF * (J - 1) + 1                 ! "FIRST" COLUMN, U

!-->     ONE-DIMENSIONAL AXIAL BAR, I_OPT = 1
    B (1, K) = GDH (1, J)                     ! DU/DX  FOR XX NORMAL
    IF ( I_OPT == 1 ) CYCLE ! TO NEXT J VALUE

!->        PLANE STRESS, PLANE STRAIN, AXISYMMETRIC, 3D
!          STRESS & STRAIN COMPONENT ORDER: XX, YY, AND XY
      L = K + 1                               ! "SECOND" COLUMN, V
      B (3, K) = GDH (2, J)                   ! DU/DY  FOR XY SHEAR
      B (2, L) = GDH (2, J)                   ! DV/DY  FOR YY NORMAL
      B (3, L) = GDH (1, J)                   ! DV/DX  FOR XY SHEAR
      IF ( I_OPT == 2 .OR. I_OPT == 3 ) CYCLE ! TO NEXT J VALUE

!->          AXISYMMETRIC ONLY
!            STRESS & STRAIN COMPONENT ORDER: XX, YY, XY, AND TT      
      IF ( I_OPT == 4 ) THEN
        IF ( R <= 0.d0 ) STOP 'R=0, BELAST'
        B (4, K) = H (J) / R                  ! U/R HOOP, ZZ NORMAL
        CYCLE ! TO NEXT J VALUE
      END IF ! AXISYMMETRIC

!->        3D SOLID ONLY
!          STRESS & STRAIN COMPONENT ORDER: XX, YY, XY, ZZ, XZ, YZ
          M = L + 1                           ! "THIRD" COLUMN, W
          B (5, K) = GDH (3, J)               ! DU/DZ  FOR XZ SHEAR
          B (6, L) = GDH (3, J)               ! DV/DZ  FOR YZ SHEAR
          B (4, M) = GDH (3, J)               ! DW/DZ  FOR ZZ NORMAL
          B (5, M) = GDH (1, J)               ! DW/DX  FOR XZ SHEAR
          B (6, M) = GDH (2, J)               ! DW/DY  FOR YZ SHEAR
  END DO
END SUBROUTINE B_MATRIX_ELASTIC
                      
SUBROUTINE B_MATRIX_PLANE_ELASTIC (NOD_PER_EL, N_SPACE,    &
                                   N_G_DOF, GDH, NS, B)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
! 2-D ELASTICITY STRAIN-DISPLACEMENT RELATIONS  (B)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
Use Precision_Module
 IMPLICIT NONE
 INTEGER,  INTENT(IN)  :: NOD_PER_EL, N_SPACE, N_G_DOF, NS
 !b xxx REAL(DP), INTENT(IN)  :: GDH (N_SPACE, NOD_PER_EL * N_G_DOF)
 REAL(DP), INTENT(IN)  :: GDH (N_SPACE, NOD_PER_EL) !b
 REAL(DP), INTENT(OUT) :: B (NS, NOD_PER_EL * N_G_DOF)
 INTEGER :: J, K, L
                               
! B          = STRAIN-DISPLACEMENT MATRIX (RETURNED)
! GDH        = GLOBAL DERIVATIVES OF ELEM INTERPOLATION FUNCTIONS
! NOD_PER_EL = NUMBER OF NODES PER ELEMENT
! N_G_DOF    = NUMBER OF PARAMETERS PER NODE
! NS         = NUMBER OF STRAINS (ROWS IN B) XX, YY, XY
! N_SPACE    = DIMENSION OF SPACE
 
!b print *,'size (GDH) ', size (GDH, 1), size (GDH, 2) !b
!b print *,'size (B  ) ', size (B  , 1), size (B  , 2) !b
  DO J = 1, NOD_PER_EL        ! ROW NUBMER
    K = N_G_DOF * (J - 1) + 1 ! FIRST COLUMN, U
    L = K + 1                 ! SECOND COLUMN, V
    B (1, K) = GDH (1, J)     ! DU/DX  FOR XX NORMAL
    B (2, K) = 0.d0
    B (3, K) = GDH (2, J)     ! DU/DY  FOR XY SHEAR
    B (1, L) = 0.d0
    B (2, L) = GDH (2, J)     ! DV/DY  FOR YY NORMAL
    B (3, L) = GDH (1, J)     ! DV/DX  FOR XY SHEAR
  END DO
END SUBROUTINE B_MATRIX_PLANE_ELASTIC

SUBROUTINE HOOKE (E, STRAIN, STRAIN_0, STRESS, N_R_B)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * *
!          STRESSES DUE TO INITIAL STRAINS
! STRESS(L) = E(L,M) * ( STRAIN(M) - STRAIN_0(M) )
! * * * * * * * * * * * * * * * * * * * * * * * * * * * *
Use Precision_Module
 IMPLICIT NONE
 REAL(DP), PARAMETER   :: ZERO = 0.d0
 INTEGER,  INTENT(IN)  :: N_R_B
 REAL(DP), INTENT(IN)  :: E (N_R_B, N_R_B)
 REAL(DP), INTENT(IN)  :: STRAIN (N_R_B), STRAIN_0 (N_R_B)
 REAL(DP), INTENT(OUT) :: STRESS (N_R_B)
                                           
! E        = CONSTITUTIVE MATRIX
! N_R_B    = NUMBER OF ROWS IN B MATRIX
! STRAIN   = MECHANICAL STRAIN VECTOR
! STRAIN_0 = INITIAL STRAIN VECTOR
! STRESS   = STRESS VECTOR
                        
  STRESS = MATMUL ( E, (STRAIN - STRAIN_0))
END SUBROUTINE HOOKE
                    
SUBROUTINE VON_MISES_STRESS (STRESS, PR, I_OPT, EFFECTIVE, NS) 
!  * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
!     VON MISES EFFECTIVE STRESS FROM STRESS VECTOR            
!  * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
Use Precision_Module
  INTEGER,  INTENT(IN)  :: NS, I_OPT
  REAL(DP), INTENT(IN)  :: STRESS (NS), PR
  REAL(DP), INTENT(OUT) :: EFFECTIVE 
  REAL(DP), PARAMETER   :: A = 0.7071067812d0, SIX = 6.d0 
  REAL(DP) :: STRESS_4

!  I_OPT     = Options: 1=1D, 2=plane stress, 3=plane strain, 
!                       4=axisymmetric stress, 5=3D solid 
!  NS        = Number of stress components                    
!  PR        = Poisson ratio                                 
!  STRESS    = Stress vector components                  
!               I_OPT = 1, Sxx, Axial
!               I_OPT = 2, Sxx  Syy  Sxy, Plane stress
!               I_OPT = 3, Sxx  Syy  Sxy  Szz, Plane strain
!               I_OPT = 4, Sxx  Syy  Sxy  Szz, Axisymmetric
!               I_OPT = 5, Sxx  Syy  Sxy  Szz  Sxz  Syz, Full 3-D
!  EFFECTIVE = Von Mises Effective Stress                  
SELECT CASE (I_OPT)

  CASE (1) !      Stress = Sxx, Axial                
    EFFECTIVE = ABS (STRESS (1)) 

  CASE (2) !      Stress = Sxx  Syy  Sxy, Plane stress 
    EFFECTIVE = A * SQRT ( (STRESS (1) - STRESS (2)) **2         &
                          + STRESS (2) **2                       &
                          + STRESS (1) **2 + SIX * STRESS (3) **2)

  CASE (3) !      Stress = Sxx  Syy  Sxy  Szz, Plane strain    
    STRESS_4  = PR * (STRESS (1) + STRESS (2)) 
    EFFECTIVE = A * SQRT ( (STRESS (1) - STRESS (2)) **2  &
                         + (STRESS (2) - STRESS_4) **2    &
                         + (STRESS_4 - STRESS (1)) **2    &
                         + SIX * STRESS (3) **2)

  CASE (4) !      Stress = Sxx  Syy  Sxy  Szz, Axisymmetric      
    EFFECTIVE = A * SQRT ( (STRESS (1) - STRESS (2)) **2 &
                         + (STRESS (2) - STRESS (4)) **2 &
                         + (STRESS (4) - STRESS (1)) **2 &
                         + SIX * STRESS (3) **2)                      

  CASE (5) !      Stress = Sxx   Syy  Sxy  Szz  Sxz  Syz, Full 3-D 
    EFFECTIVE = A * SQRT ( (STRESS (1) - STRESS (2)) **2          &
                         + (STRESS (2) - STRESS (4)) **2          &
                         + (STRESS (4) - STRESS (1)) **2          &
                         + SIX * (STRESS (3) **2 + STRESS (5) **2 &
                         + STRESS (6) **2))              
  CASE DEFAULT
    STOP 'INVALID I_OPT, VON_MISES_STRESS'
END SELECT
END SUBROUTINE VON_MISES_STRESS 

SUBROUTINE AXISYMMETRIC_E_MATRIX (Y_M, P_R, E)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
! AXISYMMETRIC CONSTITUTIVE MATRIX DEFINITION, N_R_B = 4
! STRESS & STRAIN COMPONENT ORDER: RR, ZZ, RZ, AND TT      
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants   ! for N_R_B
 IMPLICIT NONE
 REAL(DP), INTENT(IN)  :: Y_M               ! YOUNG'S MODULUS
 REAL(DP), INTENT(IN)  :: P_R               ! POISSON'S RATIO
 REAL(DP), INTENT(OUT) :: E (N_R_B, N_R_B)  ! CONSTITUTIVE
 REAL(DP) :: C1, C2, C3, S_M

! N_R_B = NUMBER OF ROWS IN B AND E MATRICES
! S_M   = SHEAR MODULUS

  IF ( N_R_B /= 4 ) STOP 'INVALID b_rows IN AXISYMMETRIC_E_MATRIX'
  S_M = 0.5d0 * Y_M / ( 1.0d0 + P_R)
  C1  = Y_M /(1.d0 + P_R)/(1.d0 - P_R - P_R)
  C2  = C1 * P_R  ; C3 = C1 * (1.d0 - P_R)

  E (:, 3) = 0.d0 ; E (3, :) = 0.d0
  E (1, 1) = C3   ; E (2, 1) = C2   ; E (4, 1) = C2  
  E (1, 2) = C2   ; E (2, 2) = C3   ; E (4, 2) = C2  
  E (1, 4) = C2   ; E (2, 4) = C2   ; E (4, 4) = C3 
  E (3, 3) = S_M 
!b if ( I_BUG > 0 ) call rprint (E, N_R_B, N_R_B, 1) !b
END SUBROUTINE AXISYMMETRIC_E_MATRIX 

SUBROUTINE AXISYMMETRIC_E_MATRIX_GEN (Y_M, P_R, S_M, E)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
! AXISYMMETRIC CONSTITUTIVE MATRIX, INDEPENDENT SHEAR MODULUS
! STRESS & STRAIN COMPONENT ORDER: RR, ZZ, RZ, AND TT      
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants   ! for N_R_B
 IMPLICIT NONE
 REAL(DP), INTENT(IN)  :: Y_M               ! YOUNG'S MODULUS
 REAL(DP), INTENT(IN)  :: P_R               ! POISSON'S RATIO
 REAL(DP), INTENT(IN)  :: S_M               ! SHEAR MODULUS
 REAL(DP), INTENT(OUT) :: E (N_R_B, N_R_B)  ! CONSTITUTIVE
 REAL(DP) :: C1, C2, C3

! N_R_B = NUMBER OF ROWS IN B AND E MATRICES

  IF ( N_R_B /= 4 ) STOP 'INVALID b_rows IN AXISYMMETRIC_E_MATRIX_GEN'
  C1 = Y_M /(1.d0 + P_R)/(1.d0 - P_R - P_R)
  C2 = C1 * P_R  ; C3 = C1 * (1.d0 - P_R)

  E (:, 3) = 0.d0 ; E (3, :) = 0.d0
  E (1, 1) = C3   ; E (2, 1) = C2   ; E (4, 1) = C2  
  E (1, 2) = C2   ; E (2, 2) = C3   ; E (4, 2) = C2  
  E (1, 4) = C2   ; E (2, 4) = C2   ; E (4, 4) = C3 
  E (3, 3) = S_M 
!b if ( I_BUG > 0 ) call rprint (E, N_R_B, N_R_B, 1) !b
END SUBROUTINE AXISYMMETRIC_E_MATRIX_GEN 

SUBROUTINE PLANE_STRESS_E_MATRIX (Y, PR, E)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
!     PLANE_STRESS CONSTITUTIVE MATRIX DEFINITION
! STRESS & STRAIN COMPONENT ORDER: XX, YY, XY, SO N_R_B = 3
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants   ! for N_R_B
 IMPLICIT NONE
 REAL(DP), INTENT(IN)  :: Y, PR
 REAL(DP), INTENT(OUT) :: E (N_R_B, N_R_B)
 REAL(DP) :: CONST, SM

! E     = CONSTITUTIVE MATRIX
! N_R_B = NUMBER OF ROWS IN B AND E MATRICES
! PR    = POISSON'S RATIO                                       
! SM    = SHEAR MODULUS
! Y     = YOUNG'S MODULUS OF ELASTICITY                         

  IF ( N_R_B /= 3 ) STOP 'INVALID b_rows IN PLANE_STRESS_E_MATRIX'

  SM    = 0.5d0 * Y / ( 1.0d0 + PR)
  CONST = Y / (1.0d0 - PR * PR) 

  E (1, 1) = CONST      ; E (2, 1) = CONST * PR ; E (3, 1) = 0.d0
  E (1, 2) = CONST * PR ; E (2, 2) = CONST      ; E (3, 2) = 0.d0
  E (1, 3) = 0.d0       ; E (2, 3) = 0.d0       ; E (3, 3) = SM 
!b if ( I_BUG > 0 ) print *, 'plane stress E ' !b
!b if ( I_BUG > 0 ) call rprint (E, N_R_B, N_R_B, 1) !b
END SUBROUTINE PLANE_STRESS_E_MATRIX 

SUBROUTINE PLANE_STRESS_E_MATRIX_GEN (Y, PR, SM, E)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
!     PLANE_STRESS CONSTITUTIVE MATRIX DEFINITION
!       WITH INDEPENDENT SHEAR MODULUS, SM
! STRESS & STRAIN COMPONENT ORDER: XX, YY, XY, SO N_R_B = 3
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants   ! for N_R_B
 IMPLICIT NONE
 REAL(DP), INTENT(IN)  :: Y, PR, SM
 REAL(DP), INTENT(OUT) :: E (N_R_B, N_R_B)
 REAL(DP) :: CONST

! E     = CONSTITUTIVE MATRIX
! N_R_B = NUMBER OF ROWS IN B AND E MATRICES
! PR    = POISSON'S RATIO                                       
! SM    = SHEAR MODULUS
! Y     = YOUNG'S MODULUS OF ELASTICITY                         

  IF ( N_R_B /= 3 ) STOP 'INVALID b_rows IN PLANE_STRESS_E_MATRIX_GEN'

  CONST = Y / (1.0d0 - PR * PR) 

  E (1, 1) = CONST      ; E (2, 1) = CONST * PR ; E (3, 1) = 0.d0
  E (1, 2) = CONST * PR ; E (2, 2) = CONST      ; E (3, 2) = 0.d0
  E (1, 3) = 0.d0       ; E (2, 3) = 0.d0       ; E (3, 3) = SM 
!b if ( I_BUG > 0 ) print *, 'plane stress E generalized ' !b
!b if ( I_BUG > 0 ) call rprint (E, N_R_B, N_R_B, 1) !b
END SUBROUTINE PLANE_STRESS_E_MATRIX_GEN 

SUBROUTINE PLANE_STRAIN_E_MATRIX (Y, PR, E)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
!     PLANE_STRAIN CONSTITUTIVE MATRIX DEFINITION
! STRESS & STRAIN COMPONENT ORDER: XX, YY, XY, SO N_R_B = 3
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants   ! for N_R_B
 IMPLICIT NONE
 REAL(DP), INTENT(IN)  :: Y, PR
 REAL(DP), INTENT(OUT) :: E (N_R_B, N_R_B)
 REAL(DP) :: C_1, C_2, SM  

! E     = CONSTITUTIVE MATRIX
! N_R_B = NUMBER OF ROWS IN B AND E MATRICES
! PR    = POISSON'S RATIO                                       
! SM    = SHEAR MODULUS
! Y     = YOUNG'S MODULUS OF ELASTICITY                         

  IF ( N_R_B /= 3 )  STOP 'INVALID b_rows IN PLANE_STRAIN_E_MATRIX'
  IF ( PR == 0.5d0 ) STOP 'INVALID PROPERTY IN PLANE_STRAIN_E_MATRIX'

  C_1 = Y * (1.d0 - PR) / ((1.d0 + PR)*(1.d0 - PR - PR))
  C_2 = Y *  PR / ((1.d0 + PR)*(1.d0 - PR - PR))
  SM  = 0.5d0 * Y / ( 1.0d0 + PR)

  E (1, 1) = C_1        ; E (2, 1) = C_2        ; E (3, 1) = 0.d0
  E (1, 2) = C_2        ; E (2, 2) = C_1        ; E (3, 2) = 0.d0
  E (1, 3) = 0.d0       ; E (2, 3) = 0.d0       ; E (3, 3) = SM 
END SUBROUTINE PLANE_STRAIN_E_MATRIX 

SUBROUTINE ELASTIC_B_AXISYMMETRIC (DGH, R, B)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
! AXISYMMETRIC ELASTICITY STRAIN-DISPLACEMENT RELATIONS  (B)
! STRESS & STRAIN COMPONENT ORDER: RR, ZZ, RZ, AND TT      
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants ! for DP, N_R_B, N_G_DOF, N_SPACE
Use Elem_Type_Data   ! for LT_FREE, LT_N, H (LT_N)
 IMPLICIT NONE
 REAL(DP), INTENT(IN)  :: DGH (N_SPACE, LT_N) ! Gradients
 REAL(DP), INTENT(IN)  :: R                   ! Radius
 REAL(DP), INTENT(OUT) :: B   (N_R_B, LT_N * N_G_DOF)
 INTEGER :: J, K, L
                               
! B       = STRAIN-DISPLACEMENT MATRIX (RETURNED)
! DGH     = GLOBAL DERIVATIVES OF H
! H       = ELEMENT INTERPOLATION FUNCTIONS
! LT_N    = NUMBER OF NODES PER ELEMENT
! N_G_DOF = NUMBER OF PARAMETERS PER NODE = 2 here
! N_R_B   = NUMBER OF STRAINS (ROWS IN B) = 4: XX, YY, XY, HOOP
! N_SPACE = DIMENSION OF SPACE = 2 here
! R,Z,T   DENOTE RADIAL, AXIAL, CIRCUMFERENCE
 
  B = 0.d0  
  DO J = 1, LT_N             ! ROW NUMBER
    K = N_G_DOF * (J - 1) + 1      ! FIRST COLUMN,  U
    L = K + 1                      ! SECOND COLUMN, V
    B (1, K) = DGH (1, J)          ! DU/DX  FOR XX NORMAL
    B (3, K) = DGH (2, J)          ! DU/DY  FOR XY SHEAR
    IF ( R <= 0.d0 ) STOP 'R=0, IN ELASTIC_B_AXISYMMETRIC'
    B (4, K) = H (J) / R           ! U/R HOOP, ZZ NORMAL
    B (2, L) = DGH (2, J)          ! DV/DY  FOR YY NORMAL
    B (3, L) = DGH (1, J)          ! DV/DX  FOR XY SHEAR
  END DO
END SUBROUTINE ELASTIC_B_AXISYMMETRIC
                           
SUBROUTINE ELASTIC_B_BAR (DGH, B)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
! ELASTIC BAR STRAIN-DISPLACEMENT RELATION  (B)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants ! for DP, N_R_B, N_G_DOF, N_SPACE
Use Elem_Type_Data   ! for LT_FREE, LT_N
 IMPLICIT NONE
 REAL(DP), INTENT(IN)  :: DGH (N_SPACE, LT_N)
 REAL(DP), INTENT(OUT) :: B   (N_R_B, LT_N * N_G_DOF)
 INTEGER :: J, K
                            
! B       = STRAIN-DISPLACEMENT MATRIX (RETURNED)
! DGH     = GLOBAL DERIVATIVES OF ELEM INTERPOLATION FUNCTIONS
! LT_N    = NUMBER OF NODES PER ELEMENT
! N_G_DOF = NUMBER OF PARAMETERS PER NODE
! N_R_B   = NUMBER OF STRAINS (ROWS IN B) = 1
! N_SPACE = DIMENSION OF SPACE
 
  DO J = 1, LT_N             ! ROW NUMBER
    K = N_G_DOF * (J - 1) + 1      ! FIRST COLUMN, U
    B (1, K) = DGH (1, J)          ! DU/DX  FOR XX NORMAL
  END DO
END SUBROUTINE ELASTIC_B_BAR

SUBROUTINE ELASTIC_B_CYLINDER (DGH, R, B)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
! ELASTIC  CYLINDER STRAIN-DISPLACEMENT RELATIONS  (B)
!      STRESS & STRAIN COMPONENT ORDER: RR, AND TT      
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants ! for DP, N_R_B, N_G_DOF, N_SPACE
Use Elem_Type_Data   ! for LT_FREE, LT_N, H (LT_N)
 IMPLICIT NONE
 REAL(DP), INTENT(IN)  :: DGH (1, LT_N), R
 REAL(DP), INTENT(OUT) :: B   (2, LT_N)
 INTEGER :: J
                               
! B    = STRAIN-DISPLACEMENT MATRIX (RETURNED)
! DGH  = GLOBAL DERIVATIVES OF H
! H    = ELEMENT INTERPOLATION FUNCTIONS
! LT_N = NUMBER OF NODES PER ELEMENT
! R,T  DENOTE RADIAL, CIRCUMFERENCE
 
  B = 0.d0  
  DO J = 1, LT_N             ! ROW NUMBER
    B (1, J) = DGH (1, J)          ! DU/DR  FOR RR NORMAL
    IF ( R <= 0.d0 ) STOP 'R=0, IN ELASTIC_B_CYLINDER'
    B (2, J) = H (J) / R           ! U/R HOOP
  END DO
END SUBROUTINE ELASTIC_B_CYLINDER

SUBROUTINE ELASTIC_B_PLANAR (DGH, B)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
! 2-D ELASTICITY STRAIN-DISPLACEMENT RELATIONS     (B)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants ! for DP, N_R_B, N_G_DOF, N_SPACE
Use Elem_Type_Data   ! for LT_FREE, LT_N
 IMPLICIT NONE
 REAL(DP), INTENT(IN)  :: DGH (N_SPACE, LT_N) 
 REAL(DP), INTENT(OUT) :: B   (N_R_B, LT_N * N_G_DOF)
 INTEGER :: J, K, L
                               
! B       = STRAIN-DISPLACEMENT MATRIX (RETURNED)
! DGH     = GLOBAL DERIVATIVES OF ELEM INTERPOLATION FUNCTIONS   
! LT_N    = NUMBER OF NODES PER ELEMENT
! N_G_DOF = NUMBER OF PARAMETERS PER NODE = 2
! N_R_B   = NUMBER OF STRAINS (ROWS IN B) = 3: XX, YY, XY
! N_SPACE = DIMENSION OF SPACE = 2 here
 
  DO J = 1, LT_N        ! ROW NUBMER
    K = N_G_DOF * (J - 1) + 1 ! FIRST COLUMN, U
    L = K + 1                 ! SECOND COLUMN, V
    B (1, K) = DGH (1, J)     ! DU/DX  FOR XX NORMAL
    B (2, K) = 0.d0
    B (3, K) = DGH (2, J)     ! DU/DY  FOR XY SHEAR
    B (1, L) = 0.d0
    B (2, L) = DGH (2, J)     ! DV/DY  FOR YY NORMAL
    B (3, L) = DGH (1, J)     ! DV/DX  FOR XY SHEAR
  END DO
END SUBROUTINE ELASTIC_B_PLANAR
                           
SUBROUTINE ELASTIC_B_SOLID (DGH, B)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
! 3-D ELASTICITY STRAIN-DISPLACEMENT RELATIONS  (B)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants ! for DP, N_R_B, N_G_DOF, N_SPACE
Use Elem_Type_Data   ! for LT_FREE, LT_N
 IMPLICIT NONE
 REAL(DP), INTENT(IN)  :: DGH (N_SPACE, LT_N)
 REAL(DP), INTENT(OUT) :: B   (N_R_B, LT_N * N_G_DOF)
 INTEGER  :: J, K, L, M
                                  
! B       = STRAIN-DISPLACEMENT MATRIX (RETURNED)
! DGH     = GLOBAL DERIVATIVES OF ELEM INTERPOLATION FUNCTIONS
! LT_N    = NUMBER OF NODES PER ELEMENT
! N_G_DOF = NUMBER OF PARAMETERS PER NODE
! N_R_B   = NUMBER OF STRAINS (ROWS IN B) XX, YY, XY, ZZ, XZ, YZ
! N_SPACE = DIMENSION OF SPACE
 
  B = 0.d0
  DO J = 1, LT_N              ! ROW NUMBER
    K = N_G_DOF * (J - 1) + 1 ! FIRST COLUMN,  U TERMS
    L = K + 1                 ! SECOND COLUMN, V TERMS
    M = L + 1                 ! THIRD COLUMN,  W TERMS
    B (1, K) = DGH (1, J)     ! DU/DX  FOR XX NORMAL
    B (3, K) = DGH (2, J)     ! DU/DY  FOR XY SHEAR
    B (5, K) = DGH (3, J)     ! DU/DZ  FOR XZ SHEAR
    B (2, L) = DGH (2, J)     ! DV/DY  FOR YY NORMAL
    B (3, L) = DGH (1, J)     ! DV/DX  FOR XY SHEAR
    B (6, L) = DGH (3, J)     ! DV/DZ  FOR YZ SHEAR
    B (4, M) = DGH (3, J)     ! DW/DZ  FOR ZZ NORMAL
    B (5, M) = DGH (1, J)     ! DW/DX  FOR XZ SHEAR
    B (6, M) = DGH (2, J)     ! DW/DY  FOR YZ SHEAR
  END DO
END SUBROUTINE ELASTIC_B_SOLID

SUBROUTINE E_PLANE_STRESS (E)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
!     PLANE_STRESS CONSTITUTIVE MATRIX DEFINITION
! STRESS & STRAIN COMPONENT ORDER: XX, YY, XY, SO N_R_B = 3
! PROPERTY ORDER: 1-YOUNG'S MODULUS, 2-POISSON'S RATIO
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants    ! for DP, N_R_B
Use Sys_Properties_Data ! for function GET_REAL_LP
 IMPLICIT NONE
 REAL(DP), INTENT(OUT) :: E (N_R_B, N_R_B)
 REAL(DP) :: CONST, P_R, S_M, Y_M

! E     = CONSTITUTIVE MATRIX
! N_R_B = NUMBER OF ROWS IN B AND E MATRICES
! P_R   = POISSON'S RATIO                                       
! S_M   = SHEAR MODULUS
! Y_M   = YOUNG'S MODULUS OF ELASTICITY                         

  IF ( N_R_B  /= 3 ) STOP 'INVALID b_rows IN E_PLANE_STRESS'
  IF ( EL_REAL < 2 ) STOP 'KEYWORD el_real < 2 IN E_PLANE_STRAIN'
  Y_M = GET_REAL_LP (1) ; P_R = GET_REAL_LP (2)

  S_M   = 0.5d0 * Y_M / (1 + P_R)
  CONST = Y_M / (1 - P_R * P_R) 

  E (1, 1) = CONST       ; E (2, 1) = CONST * P_R ; E (3, 1) = 0.d0
  E (1, 2) = CONST * P_R ; E (2, 2) = CONST       ; E (3, 2) = 0.d0
  E (1, 3) = 0.d0        ; E (2, 3) = 0.d0        ; E (3, 3) = S_M 
END SUBROUTINE E_PLANE_STRESS 

SUBROUTINE E_PLANE_STRAIN (E)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
!     PLANE_STRAIN CONSTITUTIVE MATRIX DEFINITION
! STRESS & STRAIN COMPONENT ORDER: XX, YY, XY, SO N_R_B = 3
! PROPERTY ORDER: 1-YOUNG'S MODULUS, 2-POISSON'S RATIO
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants    ! for DP, N_R_B
Use Sys_Properties_Data ! for function GET_REAL_LP
 IMPLICIT NONE
 REAL(DP), INTENT(OUT) :: E (N_R_B, N_R_B)
 REAL(DP) :: C_1, C_2, P_R, S_M, Y_M

! E     = CONSTITUTIVE MATRIX
! N_R_B = NUMBER OF ROWS IN B AND E MATRICES
! P_R   = POISSON'S RATIO                                       
! S_M   = SHEAR MODULUS
! Y_M   = YOUNG'S MODULUS OF ELASTICITY                         

  IF ( N_R_B  /= 3 )  STOP 'INVALID b_rows IN E_PLANE_STRAIN'
  IF ( EL_REAL < 2 ) STOP 'KEYWORD el_real < 2 IN E_PLANE_STRAIN'
  Y_M = GET_REAL_LP (1) ; P_R = GET_REAL_LP (2)
  IF ( P_R == 0.5d0 ) STOP "INVALID POISSON'S RATIO IN E_PLANE_STRAIN"

  C_1 = Y_M * (1 - P_R) / ((1 + P_R)*(1 - P_R - P_R))
  C_2 = Y_M *  P_R / ((1 + P_R)*(1 - P_R - P_R))
  S_M = 0.5d0 * Y_M / ( 1 + P_R)

  E (1, 1) = C_1     ; E (2, 1) = C_2     ; E (3, 1) = 0.d0
  E (1, 2) = C_2     ; E (2, 2) = C_1     ; E (3, 2) = 0.d0
  E (1, 3) = 0.d0    ; E (2, 3) = 0.d0    ; E (3, 3) = S_M 
END SUBROUTINE E_PLANE_STRAIN 

SUBROUTINE E_AXISYMMETRIC_STRESS (E)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
! AXISYMMETRIC CONSTITUTIVE MATRIX DEFINITION, N_R_B = 4
! STRESS & STRAIN COMPONENT ORDER: RR, ZZ, RZ, AND TT      
! PROPERTY ORDER: 1-YOUNG'S MODULUS, 2-POISSON'S RATIO
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants    ! for DP, N_R_B
Use Sys_Properties_Data ! for function GET_REAL_LP
 IMPLICIT NONE
 REAL(DP), INTENT(OUT) :: E (N_R_B, N_R_B)  ! CONSTITUTIVE
 REAL(DP) :: C_1, C_2, C_3, P_R, S_M, Y_M

! N_R_B = NUMBER OF ROWS IN B AND E MATRICES
! P_R   = POISSON'S RATIO
! S_M   = SHEAR MODULUS
! Y_M   = YOUNG'S MODULUS OF ELASTICITY	

  IF ( N_R_B  /= 4 ) STOP 'INVALID b_rows IN E_AXISYMMETRIC_STRESS'
  IF ( EL_REAL < 2 ) STOP 'KEYWORD el_real < 2 IN E_AXISYMMETRIC_STRESS'
  Y_M = GET_REAL_LP (1) ; P_R = GET_REAL_LP (2)
  S_M = 0.5d0 * Y_M / ( 1 + P_R)
  C_1 = Y_M /(1 + P_R)/(1 - P_R - P_R)
  C_2 = C_1 * P_R  ; C_3 = C_1 * (1 - P_R)

  E (:, 3) = 0.d0  ; E (3, :) = 0.d0
  E (1, 1) = C_3   ; E (2, 1) = C_2   ; E (4, 1) = C_2  
  E (1, 2) = C_2   ; E (2, 2) = C_3   ; E (4, 2) = C_2  
  E (1, 4) = C_2   ; E (2, 4) = C_2   ; E (4, 4) = C_3 
  E (3, 3) = S_M 
END SUBROUTINE E_AXISYMMETRIC_STRESS 

SUBROUTINE E_SOLID_STRESS (E)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
!     SOLID_STRESS CONSTITUTIVE MATRIX DEFINITION
! STRESS COMPONENT ORDER: XX, YY, XY, ZZ, XZ, YZ, SO N_R_B = 6
! PROPERTY ORDER: 1-YOUNG'S MODULUS, 2-POISSON'S RATIO
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants    ! for DP, N_R_B
Use Sys_Properties_Data ! for function GET_REAL_LP
 IMPLICIT NONE
 REAL(DP), INTENT(OUT) :: E (N_R_B, N_R_B)
 REAL(DP) :: C_1, C_2, P_R, S_M, Y_M

! E     = CONSTITUTIVE MATRIX
! N_R_B = NUMBER OF ROWS IN B AND E MATRICES
! P_R   = POISSON'S RATIO                                       
! S_M   = SHEAR MODULUS
! Y_M   = YOUNG'S MODULUS OF ELASTICITY                         

  IF ( N_R_B  /= 6 ) STOP 'INVALID b_rows IN E_SOLID_STRESS'
  IF ( EL_REAL < 2 ) STOP 'KEYWORD el_real < 2 IN E_SOLID_STRESS'
  Y_M = GET_REAL_LP (1) ; P_R = GET_REAL_LP (2)
  
  C_1  = Y_M * (1 - P_R) / (1 + P_R) / (1 - P_R - P_R)
  C_2  = C_1 * P_R / (1 - P_R) ; S_M  = 0.5d0 * Y_M / (1 + P_R)

  E = 0.d0 ! Mostly
  E (1, 1) = C_1 ; E (2, 1) = C_2 ; E (4, 1) = C_2  ! normals
  E (1, 2) = C_2 ; E (2, 2) = C_1 ; E (4, 2) = C_2  ! normals
  E (1, 4) = C_2 ; E (2, 4) = C_2 ; E (4, 4) = C_1  ! normals
  E (3, 3) = S_M ; E (5, 5) = S_M ; E (6, 6) = S_M  ! shears
END SUBROUTINE E_SOLID_STRESS

SUBROUTINE ELASTICITY_E_MATRIX (E)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
!      ELASTICITY CONSTITUTIVE MATRIX DEFINITION
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants    ! for DP, N_R_B, N_SPACE, 
                        ! AXISYMMETRIC, PLANE_STRAIN
 IMPLICIT NONE
 REAL(DP), INTENT(OUT) :: E (N_R_B, N_R_B)  ! CONSTITUTIVE
! N_R_B = NUMBER OF ROWS IN B AND E MATRICES
! E     = SYMMETRIC CONSTITUTIVE MATRIX

IF ( ANISOTROPIC ) THEN
  STOP 'ELASTICITY_E_MATRIX: ANISOTROPIC ELASTICITY E UNAVAILABLE'
! SELECT CASE (N_SPACE)
!   CASE (1) ! 1-D or cylinder
!     IF ( AXISYMMETRIC ) THEN      ! Components RR & TT
!       CALL ANISOTROPIC_CYL_STRESS (E)
!     ELSE                          ! Component XX
!       CALL E_BAR_STRESS (E)       ; END IF
!   CASE (2) ! 2-D or axisymmetric
!     IF ( AXISYMMETRIC ) THEN     ! Components RR ZZ RZ TT
!       CALL ANISOTROPIC_AXI_STRESS (E)
!     ELSE IF ( PLANE_STRAIN ) THEN ! Components XX YY XY
!       CALL ANISOTROPIC_2D_STRAIN (E)
!     ELSE ! plane stress            Components XX YY XY
!       CALL ANISOTROPIC_2D_STRESS (E)      ; END IF
!   CASE (3) ! 3-D
!     CALL ANISOTROPIC_3D_STRESS (E) ! XX YY XY ZZ ZX ZY
! END SELECT
ELSE ! Isotropic

  SELECT CASE (N_SPACE)
    CASE (1) ! 1-D or cylinder
      IF ( AXISYMMETRIC ) THEN      ! Components RR & TT
      ! CALL E_CYLINDER_STRESS (E)
      ELSE                          ! Component XX
      ! CALL E_BAR_STRESS (E)       
      END IF
    CASE (2) ! 2-D or axisymmetric
      IF ( AXISYMMETRIC ) THEN     ! Components RR ZZ RZ TT
        CALL E_AXISYMMETRIC_STRESS (E)
      ELSE IF ( PLANE_STRAIN ) THEN ! Components XX YY XY
        CALL E_PLANE_STRAIN (E)
      ELSE ! plane stress            Components XX YY XY
        CALL E_PLANE_STRESS (E)      ; END IF
    CASE (3) ! 3-D
      CALL E_SOLID_STRESS (E)        ! XX YY XY ZZ ZX ZY
  END SELECT
END IF ! anisotropic material
END SUBROUTINE ELASTICITY_E_MATRIX

SUBROUTINE ELASTICITY_B_MATRIX (DGH, XYZ, B)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
!   GENERAL ELASTICITY STRAIN-DISPLACEMENT RELATIONS  (B)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants ! for DP, N_R_B, N_G_DOF, N_SPACE,
                     ! AXISYMMETRIC
Use Elem_Type_Data   ! for LT_FREE, LT_N, H (LT_N)
 IMPLICIT NONE
 REAL(DP), INTENT(IN)  :: DGH (N_SPACE, LT_N) ! Gradients
 REAL(DP), INTENT(IN)  :: XYZ (N_SPACE)       ! Coordinates
 REAL(DP), INTENT(OUT) :: B   (N_R_B, LT_N * N_G_DOF)
                               
! B       = STRAIN-DISPLACEMENT MATRIX (RETURNED)
! DGH     = GLOBAL DERIVATIVES OF H
! H       = ELEMENT INTERPOLATION FUNCTIONS
! LT_N    = NUMBER OF NODES PER ELEMENT
! N_G_DOF = NUMBER OF PARAMETERS PER NODE
! N_R_B   = NUMBER OF STRAINS (ROWS IN B) 
! N_SPACE = DIMENSION OF SPACE 
! XYZ     = CURRENT POINT COORDINATES, RADIUS = XYZ (1)
 
  SELECT CASE (N_SPACE)
    CASE (1) ! 1-D or cylinder
      IF ( AXISYMMETRIC ) THEN      ! Components RR & TT
        CALL ELASTIC_B_CYLINDER (DGH, XYZ (1), B)    
      ELSE                          ! Component XX
        CALL ELASTIC_B_BAR (DGH, B) ; END IF
    CASE (2) ! 2-D or axisymmetric
      IF ( AXISYMMETRIC ) THEN     ! Components RR ZZ RZ TT
        CALL ELASTIC_B_AXISYMMETRIC (DGH, XYZ (1), B)
      ELSE ! plane stress or strain. Components XX YY XY
        CALL ELASTIC_B_PLANAR (DGH, B); END IF
    CASE (3) ! 3-D
      CALL ELASTIC_B_SOLID (DGH, B)  ! XX YY XY ZZ ZX ZY
  END SELECT
END SUBROUTINE ELASTICITY_B_MATRIX
