! copyright 2005,  J. E. Akin, all rights reserved.
SUBROUTINE DERIV_3_L (X, DH)
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-*
! FIND LOCAL DERIVATIVES FOR A 3 NODE LINE ELEMENT
!             IN NATURAL COORDINATES
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-*
Use Precision_Module
 IMPLICIT NONE
 REAL(DP), INTENT(IN)  :: X
 REAL(DP), INTENT(OUT) :: DH (3)

! DH = LOCAL DERIVATIVES OF SHAPE FUNCTIONS (SHAPE_3_L)
! X  = LOCAL COORDINATE OF POINT,    -1 TO +1
! LOCAL NODE COORD. ARE -1,0,+1    1----2----3
  DH (1) = X - 0.5d0
  DH (2) = - 2.d0 * X
  DH (3) = X + 0.5d0
END SUBROUTINE DERIV_3_L

SUBROUTINE DERIV2_3_L (X, D2H)
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-*
! FIND LOCAL SECOND DERIVATIVES FOR A 3 NODE LINE ELEMENT
!             IN NATURAL COORDINATES
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-*
Use Precision_Module
 IMPLICIT NONE
 REAL(DP), INTENT(IN)  :: X
 REAL(DP), INTENT(OUT) :: D2H (3)

! D2H = LOCAL DERIVATIVES OF SHAPE FUNCTIONS (SHAPE_3_L)
! X  = LOCAL COORDINATE OF POINT,    -1 TO +1
! LOCAL NODE COORD. ARE -1,0,+1.   1----2----3
  D2H (1) = 1.d0
  D2H (2) = - 2.d0
  D2H (3) = 1.d0
END SUBROUTINE DERIV2_3_L
                        
SUBROUTINE DERIV_6_T (S, T, DH)
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-*
! LOCAL DERIVATIVES FOR A SIX NODE UNIT TRIANGLE
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-*
Use Precision_Module
 IMPLICIT NONE
 REAL(DP), INTENT(IN)  :: S, T
 REAL(DP), INTENT(OUT) :: DH (2, 6)

! S,T = LOCAL COORDINATES, SEE SHAPE_6_T
! DH = LOCAL DERIVATIVES OF SHAPE FUNCTIONS
! DH(1,K) = DH(K)/DS,   DH(2,K)=DH(K)/DT
! NODAL COORDS : 1-(0,0)   2-(1,0)   3-(0,1)
!                4-(0.5,0)  5-(0.5,0.5)  6-(0,0.5)
  DH (1, 1) = - 3.d0 + 4.d0 * S + 4.d0 * T
  DH (1, 2) = - 1.d0 + 4.d0 * S
  DH (1, 3) = 0.d0
  DH (1, 4) = 4.d0 - 8.d0 * S - 4.d0 * T
  DH (1, 5) = 4.d0 * T
  DH (1, 6) = - 4.d0 * T
  DH (2, 1) = - 3.d0 + 4.d0 * S + 4.d0 * T
  DH (2, 2) = 0.d0
  DH (2, 3) = - 1.d0 + 4.d0 * T
  DH (2, 4) = - 4.d0 * S
  DH (2, 5) = 4.d0 * S
  DH (2, 6) = 4.d0 - 4.d0 * S - 8.d0 * T
END SUBROUTINE DERIV_6_T

SUBROUTINE DERIV2_6_T (S, T, D2H)
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-*
! LOCAL SECOND DERIVATIVES FOR A SIX NODE UNIT TRIANGLE
!                   H,SS H,ST H,TT
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-*
Use Precision_Module
 IMPLICIT NONE
 REAL(DP), INTENT(IN)  :: S, T
 REAL(DP), INTENT(OUT) :: D2H (3, 6)

! S,T = LOCAL COORDINATES, SEE SHAPE_6_T
! D2H = LOCAL DERIVATIVES OF SHAPE FUNCTIONS
! D2H(1,K) = D2H(K)/DS_DS,   D2H(2,K)=D2H(K)/DS_DT, ETC
! NODAL COORDS : 1-(0,0)   2-(1,0)   3-(0,1)
!                4-(0.5,0)  5-(0.5,0.5)  6-(0,0.5)

!  | V,AA |         | V,II |    A,B IMPLY LOCAL COORD
!  | V,AB | = (SDJ) | V,IJ |    I,J IMPLY GLOBAL COORD
!  | V,BB |         | V,JJ |

  D2H (1, 1) = 4.d0
  D2H (1, 2) = 4.d0
  D2H (1, 3) = 0.d0
  D2H (1, 4) = -8.d0
  D2H (1, 5) = 0.d0
  D2H (1, 6) = 0.d0
  D2H (2, 1) = 4.d0
  D2H (2, 2) = 0.d0
  D2H (2, 3) = 0.d0
  D2H (2, 4) = -4.d0
  D2H (2, 5) = 4.d0
  D2H (2, 6) = -4.d0
  D2H (3, 1) = 4.d0
  D2H (3, 2) = 0.d0
  D2H (3, 3) = 4.d0
  D2H (3, 4) = 0.d0
  D2H (3, 5) = 0.d0
  D2H (3, 6) = -8.d0
END SUBROUTINE DERIV2_6_T
                        
SUBROUTINE DERIV_8_Q (S, T, DH)
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-*
! FIND LOCAL DERIVATIVES OF SHAPE FUNCTIONS FOR AN
! EIGHT NODE PARAMETRIC QUADRILATERAL ELEMENT
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-*
Use Precision_Module
 IMPLICIT NONE
 REAL(DP), INTENT(IN)  :: S, T
 REAL(DP), INTENT(OUT) :: DH (2, 8)
 REAL(DP)              :: S_P, S_M, T_P, T_M

! S,T = LOCAL COORDINATES OF POINT, SEE SHAPE_8_Q
! DH = LOCAL DERIVATIVES OF SHAPE FUNCTIONS, H
! DH(1,J) = DH(J)/DS,  DH(2,J) = DH(J)/DT
! H = SHAPE FUNCTION ARRAY
  S_P = 1.d0 + S ; S_M = 1.d0 - S
  T_P = 1.d0 + T ; T_M = 1.d0 - T
  DH (1, 1) = - 0.25d0 * T_M * (S_M + S_M + T_M - 3.d0)
  DH (2, 1) = - 0.25d0 * S_M * (T_M + S_M + T_M - 3.d0)
  DH (1, 2) = 0.25d0 * T_M * (S_P + S_P + T_M - 3.d0)
  DH (2, 2) = - 0.25d0 * S_P * (T_M + S_P + T_M - 3.d0)
  DH (1, 3) = 0.25d0 * T_P * (S_P + S_P + T_P - 3.d0)
  DH (2, 3) = 0.25d0 * S_P * (T_P + S_P + T_P - 3.d0)
  DH (1, 4) = - 0.25d0 * T_P * (S_M + S_M + T_P - 3.d0)
  DH (2, 4) = 0.25d0 * S_M * (T_P + S_M + T_P - 3.d0)
  DH (1, 5) = - S * T_M
  DH (2, 5) = - 0.5d0 * (1.d0 - S * S)
  DH (1, 6) = 0.5d0 * (1.d0 - T * T)
  DH (2, 6) = - T * S_P
  DH (1, 7) = - S * T_P
  DH (2, 7) = 0.5d0 * (1.d0 - S * S)
  DH (1, 8) = - 0.5d0 * (1.d0 - T * T)
  DH (2, 8) = - T * S_M
END SUBROUTINE DERIV_8_Q
                        
SUBROUTINE DCHECK (DELTA, NOD_PER_EL, N_SPACE)
! * * * * * * * * * * * * * * * * * * * * * * * * * *
! CHECKING OF THE LOCAL COORDINATE DERIVATIVES OF THE
! NOD_PER_EL SHAPE FUNCTIONS AT A LOCAL POINT FOR A C0 ELEMENT
! * * * * * * * * * * * * * * * * * * * * * * * * * *
Use Precision_Module
 IMPLICIT NONE
 INTEGER,  INTENT(IN) :: NOD_PER_EL, N_SPACE
 REAL(DP), INTENT(IN) :: DELTA (N_SPACE, NOD_PER_EL)
 REAL(DP)  :: D_SUM, TOL = 1.0D-7
 INTEGER   :: I_ERR, J, N_PRT = 6
                                 
! DELTA  = LOCAL DERIVATIVES OF SHAPE FUNCTIONS
! NOD_PER_EL      = NUMBER OF SHAPE FUNCTIONS
! N_SPACE = DIMENSION OF LOCAL SPACE
  I_ERR = 0
  DO J = 1, N_SPACE
    D_SUM = SUM (DELTA (J, 1:NOD_PER_EL))
    IF (ABS (D_SUM) > TOL) THEN
      IF (I_ERR == 0) WRITE (N_PRT, *) 'DERIVATIVES INCORRECT'
      I_ERR = 1 ; WRITE (N_PRT, * ) 'ROW, D_SUM', J, D_SUM
    END IF
  END DO
  IF (I_ERR /= 0) THEN
    CALL RPRINT (DELTA, NOD_PER_EL, N_SPACE, 1)
    WRITE (N_PRT, * ) 'END OF WARNING FROM DCHECK'
  END IF
END SUBROUTINE DCHECK
                     
SUBROUTINE DERIV (PT, DLH, NOD_PER_EL, N_SPACE, L_SHAPE, &
                  N_G_DOF)
                  !b N_G_DOF, ELEM_NODES)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
!     EVALUATE C0 INTERPOLATION LOCAL DERIVATIVES
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
Use Precision_Module
 IMPLICIT NONE
 INTEGER,  INTENT(IN)  :: NOD_PER_EL, N_SPACE, L_SHAPE, N_G_DOF
 !b INTEGER,  INTENT(IN)  :: ELEM_NODES (NOD_PER_EL)
 REAL(DP), INTENT(IN)  :: PT (N_SPACE)
 REAL(DP), INTENT(OUT) :: DLH (N_SPACE, NOD_PER_EL * N_G_DOF)
                                                             
! DLH        = LOCAL DERIVATIVES AT PT
!b ! ELEM_NODES = TOPOLOGY LIST, IF VARIABLE
! L_SHAPE    = 1-LINE, 2-TRI, 3-QUAD, 4-HEX, 5-TET, 6-WEDGE
! NOD_PER_EL = NUMBER OF NODES PER ELEMENT
! N_G_DOF    = NUMBER OF DEGREES OF FREEDOM PER NODE
! N_SPACE    = NUMBER OF SPATIAL DIMENSIONS
! PT         = LOCAL COORD OF A POINT

!  BRANCH ON SHAPE, THEN NUMBER OF NODES
  IF ( L_SHAPE == 1 ) THEN

!-->     1-D ELEMENTS
    IF ( NOD_PER_EL == 2 ) CALL DERIV_2_L (PT (1), DLH)
    IF ( NOD_PER_EL == 3 ) CALL DERIV_3_L (PT(1),DLH)
    RETURN
  ELSEIF ( L_SHAPE == 2 ) THEN

!-->     TRIANGULAR 2-D ELEMENTS
    IF ( NOD_PER_EL == 3  ) CALL DERIV_3_T (PT (1), PT (2), DLH)
!   IF ( NOD_PER_EL == 4  ) CALL DERIV_4_T (PT(1),PT(2),DLH)
    IF ( NOD_PER_EL == 6  ) CALL DERIV_6_T (PT (1), PT (2), DLH)
!   IF ( NOD_PER_EL == 7  ) CALL DERIV_7_T (PT(1),PT(2),DLH)
    IF ( NOD_PER_EL == 10 ) CALL DERIV_10_T (PT(1),PT(2),DLH)
    IF ( NOD_PER_EL == 15 ) CALL DERIV_15_T (PT(1),PT(2),DLH)
    RETURN
  ELSEIF ( L_SHAPE == 3 ) THEN

!-->     QUADRILATERAL 2-D ELEMENTS
    IF ( NOD_PER_EL == 4  ) CALL DERIV_4_Q  (PT (1), PT (2), DLH)
    IF ( NOD_PER_EL == 8  ) CALL DERIV_8_Q  (PT (1), PT (2), DLH)
    IF ( NOD_PER_EL == 9  ) CALL DERIV_9_Q  (PT (1), PT (2), DLH)
    IF ( NOD_PER_EL == 12 ) CALL DERIV_12_Q (PT (1), PT (2), DLH)
    IF ( NOD_PER_EL == 16 ) CALL DERIV_16_Q (PT (1), PT (2), DLH)
    IF ( NOD_PER_EL == 17 ) CALL DERIV_17_Q (PT (1), PT (2), DLH)
!   IF ( NOD_PER_EL == 25 ) CALL DERIV_25_Q (PT (1), PT (2), DLH)
    RETURN
  ELSEIF ( L_SHAPE == 4) THEN

!-->     HEXAHEDRA 3-D ELEMENTS
    IF ( NOD_PER_EL == 8  ) CALL DERIV_8_H  (PT (1), PT (2), PT (3), DLH)
    IF ( NOD_PER_EL == 20 ) CALL DERIV_20_H (PT (1), PT (2), PT (3), DLH)
!b !   IF ( NOD_PER_EL == 20 ) CALL DERIV_8_20_H (PT(1),PT(2),PT(3),&
!b !                                DLH,ELEM_NODES)
!   IF ( NOD_PER_EL == 27 ) CALL DERIV_27_H (PT(1),PT(2),PT(3),DLH)
!   IF ( NOD_PER_EL == 32 ) CALL DERIV_32_H (PT(1),PT(2),PT(3),DLH)
    RETURN
  ELSEIF ( L_SHAPE == 5) THEN

!-->     TETRAHEDRA 3-D ELEMENTS (PYRAMIDS)
    IF ( NOD_PER_EL == 4  ) CALL DERIV_4_P  (PT(1),PT(2),PT(3),DLH)
!   IF ( NOD_PER_EL == 10 ) CALL DERIV_10_P (PT(1),PT(2),PT(3),DLH)
!   IF ( NOD_PER_EL == 20 ) CALL DERIV_20_P (PT(1),PT(2),PT(3),DLH)
    RETURN
  ELSEIF ( L_SHAPE == 6) THEN

!-->     WEDGE 3-D ELEMENTS
    STOP 'NO WEDGE IN SHAPE'
!   IF ( NOD_PER_EL == 6  ) CALL DERIV_6_W  (PT(1),PT(2),PT(3),DLH)
!   IF ( NOD_PER_EL == 15 ) CALL DERIV_15_W (PT(1),PT(2),PT(3),DLH)
!   RETURN
  ELSEIF ( L_SHAPE == 7 ) THEN

!-->     USER SUPPLIED ELEMENT
!b !   CALL  DERIV_USER (PT(1),PT(2),PT(3),DLH,ELEM_NODES)
    STOP 'NO USER ELEMENT IN DERIV'
  ELSEIF ( L_SHAPE < 1 .OR. L_SHAPE > 7 ) THEN

!-->     UNSUPPORTED OPTION
    STOP 'UNSUPPORTED ELEMENT IN DERIV'
  END IF
END SUBROUTINE DERIV
                    
SUBROUTINE DERIV_6_T_BAR (R, S, T, DBH)
! * * * * * * * * * * * * * * * * * * * * * * * * * *
! BARACENTRIC DERIVATIVES FOR A T6 ELEMENT
! * * * * * * * * * * * * * * * * * * * * * * * * * *
Use Precision_Module
 IMPLICIT NONE
 REAL(DP), INTENT(IN)  :: R, S, T
 REAL(DP), INTENT(OUT) :: DBH (3, 6)
 REAL(DP), PARAMETER   :: ZERO = 0.D0, ONE = 1.D0, FOUR = 4.D0
                                                              
! DBH   = BARACENTRIC DERIVATIVES OF H
! R,S,T = BARACENTRIC COORDINATES, R+S+T=1
!  NODE  R  S  T       NODE   R   S   T     ORDER:
!  1     1  0  0       4     1/2 1/2  0     3
!  2     0  1  0       5      0  1/2 1/2    6 5
!  3     0  0  1       6     1/2  0  1/2    1 4 2
                                                 
!  BARACENTRIC DERIVATIVES
  DBH (1, 1) = FOUR * R - ONE
  DBH (2, 1) = ZERO
  DBH (3, 1) = ZERO
  DBH (1, 2) = ZERO
  DBH (2, 2) = FOUR * S - ONE
  DBH (3, 2) = ZERO
  DBH (1, 3) = ZERO
  DBH (2, 3) = ZERO
  DBH (3, 3) = FOUR * T - ONE
  DBH (1, 4) = ZERO
  DBH (2, 4) = FOUR * T
  DBH (3, 4) = FOUR * S
  DBH (1, 5) = FOUR * T
  DBH (2, 5) = ZERO
  DBH (3, 5) = FOUR * R
  DBH (1, 6) = FOUR * S
  DBH (2, 6) = FOUR * R
  DBH (3, 6) = ZERO
END SUBROUTINE DERIV_6_T_BAR
                            
SUBROUTINE DERIV_USER (PT, DH, NOD_PER_EL, N_SPACE, N_G_DOF)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
!  USER DEFINED ELEMENT C0 INTERPOLATION FUNCTION DERIVATIVES
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
Use Precision_Module
 IMPLICIT NONE
 INTEGER,  INTENT(IN)  :: NOD_PER_EL, N_SPACE, N_G_DOF
 REAL(DP), INTENT(IN)  :: PT (N_SPACE)
 REAL(DP), INTENT(OUT) :: DH (N_SPACE, NOD_PER_EL * N_G_DOF)
                                                            
! DH         = LOCAL DERIVATIVES OF INTERPOLATION FUNCTIONS AT PT
! NOD_PER_EL = NUMBER OF NODES PER ELEMENT
! N_G_DOF    = NUMBER OF DEGREES OF FREEDOM PER NODE
! N_SPACE    = NO OF SPATIAL DIMENSIONS
! PT         = LOCAL COORD OF A POINT
!
    DH = 0.d0
  PRINT *, 'DERIV_USER: USER ERROR, DERIVATIVES NOT SUPPLIED FOR'
  PRINT *, 'ELEMENT WITH ', NOD_PER_EL, ' NODES,'
  PRINT *, 'IN ', N_SPACE, ' DIMENSIONS, AND WITH'
  PRINT *, N_G_DOF, ' DEGREES OF FREEDOM PER NODE.'
  PRINT *, 'WARNING: CHECK VALUE OF shape KEYWORD'
  STOP 'No user defined element interpolations, DERIV_USER'
END SUBROUTINE DERIV_USER
                         
SUBROUTINE DERIV2_USER (PT, D2H, NOD_PER_EL, N_SPACE, N_G_DOF, N_2_DER)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
!  USER DEFINED ELEMENT C0 INTERPOLATION FUNCTION SECOND
!                      DERIVATIVES
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
Use Precision_Module
 IMPLICIT NONE
 INTEGER,  INTENT(IN)  :: NOD_PER_EL, N_SPACE, N_G_DOF, N_2_DER
 REAL(DP), INTENT(IN)  :: PT (N_SPACE)
 REAL(DP), INTENT(OUT) :: D2H (N_2_DER, NOD_PER_EL * N_G_DOF)
 INTEGER, SAVE :: WARN = 0
                                                            
! D2H        = LOCAL SECOND DERIVATIVES OF FUNCTIONS AT PT
! NOD_PER_EL = NUMBER OF NODES PER ELEMENT
! N_G_DOF    = NUMBER OF DEGREES OF FREEDOM PER NODE
! N_SPACE    = NO OF SPATIAL DIMENSIONS
! N_2_DER    = NUMBER OF SECOND DERIVATIVES
! PT         = LOCAL COORD OF A POINT
!
  D2H = 0.d0
  IF ( WARN == 0 ) THEN
    WARN = 1
    PRINT *, 'DERIV2_USER: USER ERROR, SECOND DERIVATIVES NOT SUPPLIED'
    PRINT *, 'FOR ELEMENT TYPE WITH ', NOD_PER_EL, ' NODES, IN'
    PRINT *, N_SPACE, ' DIMENSIONS, AND WITH ', N_G_DOF, ' DEGREES OF'
    PRINT *, ' FREEDOM PER NODE. CHECK VALUE OF shape KEYWORD'
    PRINT *, 'SECOND DERIVATIVES SET TO ZERO'
  END IF
  !b  STOP 'No user defined element interpolations, DERIV2_USER'
END SUBROUTINE DERIV2_USER
                         
SUBROUTINE DERIV_C_N_USER (PT, LENGTH, DH, NOD_PER_EL, &
                          N_SPACE, N_G_DOF)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
!  USER DEFINED ELEMENT C_N INTERPOLATION FUNCTION DERIVATIVES
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
Use Precision_Module
 IMPLICIT NONE
 INTEGER,  INTENT(IN)  :: NOD_PER_EL, N_SPACE, N_G_DOF
 REAL(DP), INTENT(IN)  :: PT (N_SPACE), LENGTH
 REAL(DP), INTENT(OUT) :: DH (N_SPACE, NOD_PER_EL * N_G_DOF)
                                                            
! DH         = LOCAL DERIVATIVES OF C1 INTERPOLATIONS AT PT
! LENGTH     = PHYSICAL LENGTH
! NOD_PER_EL = NUMBER OF NODES PER ELEMENT
! N_G_DOF    = NUMBER OF DEGREES OF FREEDOM PER NODE
! N_SPACE    = NO OF SPATIAL DIMENSIONS
! PT         = LOCAL COORD OF A POINT
!
    DH = 0.d0
    STOP 'No user defined element interpolations, DERIV_C_N_USER'
END SUBROUTINE DERIV_C_N_USER
                         
SUBROUTINE SHAPE_USER (PT, H, NOD_PER_EL, N_SPACE, N_G_DOF)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
!     USER DEFINED ELEMENT C0 INTERPOLATION FUNCTIONS
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
Use Precision_Module
 IMPLICIT NONE
 INTEGER,  INTENT(IN)  :: NOD_PER_EL, N_SPACE, N_G_DOF
 REAL(DP), INTENT(IN)  :: PT (N_SPACE)
 REAL(DP), INTENT(OUT) :: H (NOD_PER_EL * N_G_DOF)
                                                  
! H          = ELEMENT C0 INTERPOLATION FUNCTIONS AT PT
! NOD_PER_EL = NUMBER OF NODES PER ELEMENT
! N_G_DOF    = NUMBER OF DEGREES OF FREEDOM PER NODE
! N_SPACE    = NO OF SPATIAL DIMENSIONS
! PT         = LOCAL COORD OF A POINT

  H = 0.d0
  PRINT *, 'USER ERROR, INTERPOLATIONS NOT SUPPLIED FOR'
  PRINT *, 'ELEMENT WITH ', NOD_PER_EL, ' NODES,'
  PRINT *, 'IN ', N_SPACE, ' DIMENSIONS, AND WITH'
  PRINT *, N_G_DOF, ' DEGREES OF FREEDOM PER NODE.'
  PRINT *, 'WARNING: CHECK VALUE OF shape KEYWORD'
  STOP 'No user defined element interpolations, SHAPE_USER'
END SUBROUTINE SHAPE_USER
                         
SUBROUTINE SHAPE_C_N_USER (PT, LENGTH, H, NOD_PER_EL, &
                          N_SPACE, N_G_DOF)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
!  USER DEFINED ELEMENT C_N INTERPOLATION FUNCTIONS
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
Use Precision_Module
 IMPLICIT NONE
 INTEGER,  INTENT(IN)  :: NOD_PER_EL, N_SPACE, N_G_DOF
 REAL(DP), INTENT(IN)  :: PT (N_SPACE), LENGTH
 REAL(DP), INTENT(OUT) :: H (NOD_PER_EL * N_G_DOF)
                                                            
! H          = LOCAL C1 INTERPOLATION FUNCTIONS AT PT
! LENGTH     = PHYSICAL LENGTH
! NOD_PER_EL = NUMBER OF NODES PER ELEMENT
! N_G_DOF    = NUMBER OF DEGREES OF FREEDOM PER NODE
! N_SPACE    = NO OF SPATIAL DIMENSIONS
! PT         = LOCAL COORD OF A POINT

    H = 0.d0
    STOP 'No user defined element interpolations, SHAPE_C_N_USER'
END SUBROUTINE SHAPE_C_N_USER
                         
SUBROUTINE SCALAR_DERIVS (PT_L, DH_L)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
! EVALUATE C0 ELEMENT DERIVATIVE FUNCTIONS IN 1D, 2D, 3D
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants
Use Elem_Type_Data
 IMPLICIT NONE
 REAL(DP), INTENT(IN)  :: PT_L (LT_PARM)
 REAL(DP), INTENT(OUT) :: DH_L (LT_PARM, LT_N)
                                                            
! DH_L      = ELEMENT INTERPOLATION DERIVATIVES AT PT_L
! LT_N    = NUMBER OF NODES PER ELEMENT
! N_G_DOF = NUMBER OF DEGREES OF FREEDOM PER NODE
! LT_PARM = NO OF SPATIAL DIMENSIONS
! PT_L      = LOCAL COORD OF A POINT

!  BRANCH ON SPACE, THEN NUMBER OF NODES 
  SELECT CASE ( LT_PARM )
    CASE (0) ! point
      PRINT *,'WARNING, NO POINT ELEMENT DERIVATIVES'
      N_WARN = N_WARN + 1 ; 
    CASE (1) ! 1-D space

!-->     1-D ELEMENTS
      SELECT CASE ( LT_N )
        CASE (2) 
          SELECT CASE ( N_G_DOF )
            CASE (1) ; CALL DERIV_2_L  (PT_L (1), DH_L(1,:))
            CASE DEFAULT 
              CALL DERIV_USER (PT_L, DH_L, LT_N, LT_PARM, N_G_DOF)
          END SELECT
        CASE (3) ; CALL DERIV_3_L (PT_L (1), DH_L(1,:))
        CASE (4) ; CALL DERIV_4_L (PT_L (1), DH_L(1,:))
        CASE DEFAULT 
          CALL DERIV_USER (PT_L, DH_L, LT_N, LT_PARM, N_G_DOF)
      END SELECT ! number of element nodes
    CASE (2) ! 2-D space
      SELECT CASE ( LT_N )

!-->       TRIANGULAR 2-D ELEMENTS
        CASE ( 3) ; CALL DERIV_3_T  (PT_L (1), PT_L (2), DH_L)
!       CASE ( 4) ; CALL DERIV_4_T  (PT_L (1), PT_L (2), DH_L)
        CASE ( 6) ; CALL DERIV_6_T  (PT_L (1), PT_L (2), DH_L)
!       CASE ( 7) ; CALL DERIV_7_T  (PT_L (1), PT_L (2), DH_L)
        CASE (10) ; CALL DERIV_10_T (PT_L (1), PT_L (2), DH_L)
        CASE (15) ; CALL DERIV_15_T (PT_L (1), PT_L (2), DH_L)

!-->       QUADRILATERAL 2-D ELEMENTS
        CASE ( 4) ; CALL DERIV_4_Q  (PT_L (1), PT_L (2), DH_L)
        CASE ( 8) ; CALL DERIV_8_Q  (PT_L (1), PT_L (2), DH_L)
        CASE ( 9) ; CALL DERIV_9_Q  (PT_L (1), PT_L (2), DH_L)
        CASE (12) ; CALL DERIV_12_Q (PT_L (1), PT_L (2), DH_L)
        CASE (16) ; CALL DERIV_16_Q (PT_L (1), PT_L (2), DH_L)
        CASE (17) ; CALL DERIV_17_Q (PT_L (1), PT_L (2), DH_L)
!       CASE (25) ; CALL DERIV_25_Q (PT_L (1), PT_L (2), DH_L)
        CASE DEFAULT 
          CALL DERIV_USER (PT_L, DH_L, LT_N, LT_PARM, N_G_DOF)
      END SELECT ! number of element nodes
    CASE (3) ! 3-D space
      SELECT CASE ( LT_N )

!-->       HEXAHEDRA 3-D ELEMENTS
        CASE ( 8) ; CALL DERIV_8_H  (PT_L (1), PT_L (2), PT_L (3), DH_L)
        CASE (20) ; CALL DERIV_20_H (PT_L (1), PT_L (2), PT_L (3), DH_L)
!       CASE (27) ; CALL DERIV_27_H (PT_L (1), PT_L (2), PT_L (3), DH_L)
!       CASE (32) ; CALL DERIV_32_H (PT_L (1), PT_L (2), PT_L (3), DH_L)

!-->       TETRAHEDRA 3-D ELEMENTS (PYRAMIDS)
        CASE ( 4) ; CALL DERIV_4_P  (PT_L (1), PT_L (2), PT_L (3), DH_L)
!       CASE (10) ; CALL DERIV_10_P (PT_L (1), PT_L (2), PT_L (3), DH_L)
!       CASE (20) ; CALL DERIV_20_P (PT_L (1), PT_L (2), PT_L (3), DH_L)

!-->       WEDGE 3-D ELEMENTS
!       CASE ( 6) ; CALL DERIV_6_W  (PT_L (1), PT_L (2), PT_L (3), DH_L)
!       CASE (15) ; CALL DERIV_15_W (PT_L (1), PT_L (2), PT_L (3), DH_L)

!-->       USER SUPPLIED ELEMENT
        CASE DEFAULT 
          CALL DERIV_USER (PT_L, DH_L, LT_N, LT_PARM, N_G_DOF)
      END SELECT ! number of element nodes
    CASE DEFAULT
      PRINT *, LT_N, LT_PARM, N_G_DOF
      STOP 'INVALID SPACE SIZE IN SCALAR_DERIVS'
END SELECT ! number of spatial dimensions
END SUBROUTINE SCALAR_DERIVS
                         
SUBROUTINE SCALAR_2ND_DERIVS (PT_L, D2H_L)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
! EVALUATE C0 ELEMENT 2ND DERIVATIVE FUNCTIONS IN 1D, 2D, 3D
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants
Use Elem_Type_Data
 IMPLICIT NONE
 REAL(DP), INTENT(IN)  :: PT_L (LT_PARM)
 REAL(DP), INTENT(OUT) :: D2H_L (N_2_DER, LT_N)
                                                            
! D2H_L   = ELEMENT INTERPOLATION SECOND DERIVATIVES AT PT_L
! LT_N    = NUMBER OF NODES PER ELEMENT
! N_G_DOF = NUMBER OF DEGREES OF FREEDOM PER NODE
! LT_PARM = NO OF SPATIAL DIMENSIONS
! PT_L    = LOCAL COORD OF A POINT

!  BRANCH ON SPACE, THEN NUMBER OF NODES 
  SELECT CASE ( LT_PARM )
    CASE (0) ! point
      PRINT *,'WARNING, NO POINT ELEMENT 2ND DERIVATIVES'
      N_WARN = N_WARN + 1 ; 
    CASE (1) ! 1-D space

!-->     1-D ELEMENTS
!         | V,AA | = (P_AJ) | V,II |    A,B IMPLY LOCAL COORD
!                                       I,J IMPLY GLOBAL COORD
      SELECT CASE ( LT_N )
        CASE (2) 
          SELECT CASE ( N_G_DOF )
            CASE (1) ; D2H_L = 0.d0 
            CASE DEFAULT 
              CALL DERIV2_USER (PT_L, D2H_L, LT_N, LT_PARM, &
                                N_G_DOF, N_2_DER)
          END SELECT
        CASE (3) ; D2H_L (1, :) = (/ 1.d0, -2.d0, 1.d0 /)
        CASE (4) ; CALL DERIV2_4_L (PT_L (1), D2H_L(1,:))
        CASE DEFAULT 
          CALL DERIV2_USER (PT_L, D2H_L, LT_N, LT_PARM, &
                            N_G_DOF, N_2_DER)
      END SELECT ! number of element nodes
    CASE (2) ! 2-D space
!         | V,AA |          | V,II |    A,B IMPLY LOCAL COORD
!         | V,AB | = (P_AJ) | V,IJ |    I,J IMPLY GLOBAL COORD
!         | V,BB |          | V,JJ |

      SELECT CASE ( LT_N ) ! NUMBER OF NODES
!-->       TRIANGULAR 2-D ELEMENTS
        CASE ( 3) ; D2H_L = 0.d0 
!       CASE ( 4) ; CALL DERIV2_4_T  (PT_L (1), PT_L (2), D2H_L)
        CASE ( 6) ; CALL DERIV2_6_T  (PT_L (1), PT_L (2), D2H_L)
!       CASE ( 7) ; CALL DERIV2_7_T  (PT_L (1), PT_L (2), D2H_L)
!       CASE (10) ; CALL DERIV2_10_T (PT_L (1), PT_L (2), D2H_L)
        !b CASE (15) ; CALL DERIV2_15_T (PT_L (1), PT_L (2), D2H_L)

!-->       QUADRILATERAL 2-D ELEMENTS
        CASE ( 4) ; CALL DERIV2_4_Q  (PT_L (1), PT_L (2), D2H_L)
        !b CASE ( 8) ; CALL DERIV2_8_Q  (PT_L (1), PT_L (2), D2H_L)
        CASE ( 9) ; CALL DERIV2_9_Q  (PT_L (1), PT_L (2), D2H_L)
!       CASE (12) ; CALL DERIV2_12_Q (PT_L (1), PT_L (2), D2H_L)
!       CASE (16) ; CALL DERIV2_16_Q (PT_L (1), PT_L (2), D2H_L)
!       CASE (17) ; CALL DERIV2_17_Q (PT_L (1), PT_L (2), D2H_L)
!       CASE (25) ; CALL DERIV2_25_Q (PT_L (1), PT_L (2), D2H_L)
        CASE DEFAULT 
          CALL DERIV2_USER (PT_L,D2H_L,LT_N,LT_PARM,N_G_DOF,N_2_DER)
      END SELECT ! number of element nodes
    CASE (3) ! 3-D space
      SELECT CASE ( LT_N )

!-->       HEXAHEDRA 3-D ELEMENTS
        !b CASE ( 8) ; CALL DERIV2_8_H  (PT_L (1), PT_L (2), PT_L (3), &
                                      !b D2H_L)
!b       CASE (20) ; CALL DERIV2_20_H (PT_L (1), PT_L (2), PT_L (3), &
!b                                     D2H_L)
!       CASE (27) ; CALL DERIV2_27_H (PT_L (1), PT_L (2), PT_L (3), &
!                                     D2H_L)
!       CASE (32) ; CALL DERIV2_32_H (PT_L (1), PT_L (2), PT_L (3), &
!                                     D2H_L)

!-->       TETRAHEDRA 3-D ELEMENTS (PYRAMIDS)
!       CASE ( 4) ; D2H_L = 0.d0 
!                                     D2H_L)
!       CASE (10) ; CALL DERIV2_10_P (PT_L (1), PT_L (2), PT_L (3), &
!                                     D2H_L)
!       CASE (20) ; CALL DERIV2_20_P (PT_L (1), PT_L (2), PT_L (3), &
!                                     D2H_L)

!-->       WEDGE 3-D ELEMENTS
!       CASE ( 6) ; CALL DERIV2_6_W  (PT_L (1), PT_L (2), PT_L (3), &
!                                     D2H_L)
!       CASE (15)  CALL DERIV2_15_W (PT_L (1), PT_L (2), PT_L (3), &
!                                     D2H_L)

!-->       USER SUPPLIED ELEMENT
        CASE DEFAULT 
          CALL DERIV2_USER (PT_L, D2H_L, LT_N, LT_PARM, N_G_DOF, N_2_DER)
      END SELECT ! number of element nodes
    CASE DEFAULT
      PRINT *, LT_N, LT_PARM, N_G_DOF
      STOP 'INVALID SPACE SIZE IN SCALAR_2ND_DERIVS'
END SELECT ! number of spatial dimensions
END SUBROUTINE SCALAR_2ND_DERIVS
                         
SUBROUTINE ELEM_C_N_DERIVS (PT_L, LENGTH, DH_L)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
! EVALUATE C_N ELEMENT DERIVATIVE FUNCTIONS IN 1D, 2D, 3D
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants
Use Elem_Type_Data
 IMPLICIT NONE
 REAL(DP), INTENT(IN)  :: PT_L (LT_PARM), LENGTH
 REAL(DP), INTENT(OUT) :: DH_L (LT_PARM, LT_N)
                                                            
! DH_L      = ELEMENT INTERPOLATION DERIVATIVES AT PT_L
! LENGTH  = PHYSICAL LENGTH OF ELEMENT
! LT_N    = NUMBER OF NODES PER ELEMENT
! N_G_DOF = NUMBER OF DEGREES OF FREEDOM PER NODE
! LT_PARM = NO OF SPATIAL DIMENSIONS
! PT_L      = LOCAL COORD OF A POINT

!  BRANCH ON SPACE, THEN NUMBER OF NODES 
  SELECT CASE ( LT_PARM )
    CASE (0) ! point
      PRINT *,'WARNING, NO POINT ELEMENT C_N DERIVATIVES'
      N_WARN = N_WARN + 1 ; 
    CASE (1) ! 1-D space
!-->     1-D ELEMENTS
      SELECT CASE ( LT_N )
        CASE (2) ! 2 node line element
          SELECT CASE ( N_G_DOF )
            CASE (2) ; CALL DERIV_C1_L (PT_L (1), LENGTH, DH_L(1,:))
            CASE (3) ; CALL DERIV_C2_L (PT_L (1), LENGTH, DH_L(1,:))
            CASE (4) ; CALL DERIV_C3_L (PT_L (1), LENGTH, DH_L(1,:))
            CASE DEFAULT 
              CALL DERIV_C_N_USER (PT_L, LENGTH, DH_L, LT_N, &
                                  LT_PARM, N_G_DOF)
          END SELECT
        CASE DEFAULT 
          CALL DERIV_C_N_USER (PT_L, LENGTH, DH_L, LT_N,   &
                              LT_PARM, N_G_DOF)
      END SELECT ! number of element nodes
    CASE DEFAULT
      PRINT *, LT_N, LT_PARM, N_G_DOF
      STOP 'INVALID SPACE SIZE IN ELEM_C_N_DERIVS'
END SELECT ! number of spatial dimensions
END SUBROUTINE ELEM_C_N_DERIVS
                         
SUBROUTINE GEN_ELEM_DERIV (PT, DH_L, NOD_PER_EL, N_SPACE, N_G_DOF)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
! EVALUATE C0 ELEMENT DERIVATIVE FUNCTIONS IN 1D, 2D, 3D
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
Use Precision_Module
 IMPLICIT NONE
 INTEGER,  INTENT(IN)  :: NOD_PER_EL, N_SPACE, N_G_DOF
 REAL(DP), INTENT(IN)  :: PT (N_SPACE)
 REAL(DP), INTENT(OUT) :: DH_L (N_SPACE, NOD_PER_EL * N_G_DOF)
                                                            
! DH_L         = ELEMENT INTERPOLATION DERIVATIVES AT PT
! NOD_PER_EL = NUMBER OF NODES PER ELEMENT
! N_G_DOF    = NUMBER OF DEGREES OF FREEDOM PER NODE
! N_SPACE    = NO OF SPATIAL DIMENSIONS
! PT         = LOCAL COORD OF A POINT

!  BRANCH ON SPACE, THEN NUMBER OF NODES (and then N_G_DOF)
  SELECT CASE ( N_SPACE )
    CASE (1) ! 1-D space

!-->     1-D ELEMENTS
      SELECT CASE ( NOD_PER_EL )
        CASE (2) ; CALL DERIV_2_L (PT (1), DH_L(1,:))
        CASE (3) ; CALL DERIV_3_L (PT (1), DH_L(1,:))
        CASE (4) ; CALL DERIV_4_L (PT (1), DH_L(1,:))
        CASE DEFAULT 
          CALL DERIV_USER (PT, DH_L, NOD_PER_EL, N_SPACE, N_G_DOF)
      END SELECT ! number of element nodes
    CASE (2) ! 2-D space
      SELECT CASE ( NOD_PER_EL )

!-->       TRIANGULAR 2-D ELEMENTS
        CASE ( 3) ; CALL DERIV_3_T  (PT (1), PT (2), DH_L)
!       CASE ( 4) ; CALL DERIV_4_T  (PT (1), PT (2), DH_L)
        CASE ( 6) ; CALL DERIV_6_T  (PT (1), PT (2), DH_L)
!       CASE ( 7) ; CALL DERIV_7_T  (PT (1), PT (2), DH_L)
        CASE (10) ; CALL DERIV_10_T (PT (1), PT (2), DH_L)
        CASE (15) ; CALL DERIV_15_T (PT (1), PT (2), DH_L)

!-->       QUADRILATERAL 2-D ELEMENTS
        CASE ( 4) ; CALL DERIV_4_Q  (PT (1), PT (2), DH_L)
        CASE ( 8) ; CALL DERIV_8_Q  (PT (1), PT (2), DH_L)
        CASE ( 9) ; CALL DERIV_9_Q  (PT (1), PT (2), DH_L)
        CASE (12) ; CALL DERIV_12_Q (PT (1), PT (2), DH_L)
        CASE (16) ; CALL DERIV_16_Q (PT (1), PT (2), DH_L)
        CASE (17) ; CALL DERIV_17_Q (PT (1), PT (2), DH_L)
!       CASE (25) ; CALL DERIV_25_Q (PT (1), PT (2), DH_L)
        CASE DEFAULT 
          CALL DERIV_USER (PT, DH_L, NOD_PER_EL, N_SPACE, N_G_DOF)
      END SELECT ! number of element nodes
    CASE (3) ! 3-D space
      SELECT CASE ( NOD_PER_EL )

!-->       HEXAHEDRA 3-D ELEMENTS
        CASE ( 8) ; CALL DERIV_8_H  (PT (1), PT (2), PT (3), DH_L)
        CASE (20) ; CALL DERIV_20_H (PT (1), PT (2), PT (3), DH_L)
!       CASE (27) ; CALL DERIV_27_H (PT (1), PT (2), PT (3), DH_L)
!       CASE (32) ; CALL DERIV_32_H (PT (1), PT (2), PT (3), DH_L)

!-->       TETRAHEDRA 3-D ELEMENTS (PYRAMIDS)
        CASE ( 4) ; CALL DERIV_4_P  (PT (1), PT (2), PT (3), DH_L)
!       CASE (10) ; CALL DERIV_10_P (PT (1), PT (2), PT (3), DH_L)
!       CASE (20) ; CALL DERIV_20_P (PT (1), PT (2), PT (3), DH_L)

!-->       WEDGE 3-D ELEMENTS
!       CASE ( 6)  ; CALL DERIV_6_W  (PT (1), PT (2), PT (3), DH_L)
!       CASE (15) ; CALL DERIV_15_W (PT (1), PT (2), PT (3), DH_L)

!-->       USER SUPPLIED ELEMENT
        CASE DEFAULT 
          CALL DERIV_USER (PT, DH_L, NOD_PER_EL, N_SPACE, N_G_DOF)
      END SELECT ! number of element nodes
    CASE DEFAULT
      STOP 'INVALID SPACE SIZE IN GEN_ELEM_DERIV'
END SELECT ! number of spatial dimensions
END SUBROUTINE GEN_ELEM_DERIV
                         
SUBROUTINE SCALAR_SHAPES (PT_L, H_L)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
! EVALUATE C0 ELEMENT INTERPOLATION FUNCTIONS IN 1D, 2D, 3D
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants
Use Elem_Type_Data
 IMPLICIT NONE
 REAL(DP), INTENT(IN)  :: PT_L (LT_PARM)
 REAL(DP), INTENT(OUT) :: H_L  (LT_N)
                                                  
! H_L       = ELEMENT INTERPOLATION FUNCTIONS AT PT_L
! LT_N    = NUMBER OF NODES PER ELEMENT
! N_G_DOF = NUMBER OF DEGREES OF FREEDOM PER NODE
! LT_PARM = NO OF SPATIAL DIMENSIONS
! PT_L      = LOCAL COORD OF A POINT

!  BRANCH ON SPACE, THEN NUMBER OF NODES 
  SELECT CASE ( LT_PARM )
    CASE (0) ! point
      PRINT *,'WARNING, NO POINT ELEMENT SHAPES'     
      N_WARN = N_WARN + 1 ; H = 1.d0
    CASE (1) ! 1-D space

!-->     1-D ELEMENTS
      SELECT CASE ( LT_N )
        CASE (2) 
          SELECT CASE ( N_G_DOF )
            CASE (1) ; CALL SHAPE_2_L  (PT_L (1), H_L)
            CASE DEFAULT 
              CALL SHAPE_USER (PT_L, H_L, LT_N, LT_PARM, N_G_DOF)
          END SELECT
        CASE (3) ; CALL SHAPE_3_L (PT_L (1), H_L)
        CASE (4) ; CALL SHAPE_4_L (PT_L (1), H_L)
        CASE DEFAULT 
          CALL SHAPE_USER (PT_L, H_L, LT_N, LT_PARM, N_G_DOF)
      END SELECT ! number of element nodes
    CASE (2) ! 2-D space
      SELECT CASE ( LT_N )

!-->       TRIANGULAR 2-D ELEMENTS
        CASE ( 3) ; CALL SHAPE_3_T  (PT_L (1), PT_L (2), H_L)
!       CASE ( 4) ; CALL SHAPE_4_T  (PT_L (1), PT_L (2), H_L)
        CASE ( 6) ; CALL SHAPE_6_T  (PT_L (1), PT_L (2), H_L)
!       CASE ( 7) ; CALL SHAPE_7_T  (PT_L (1), PT_L (2), H_L)
        CASE (10) ; CALL SHAPE_10_T (PT_L (1), PT_L (2), H_L)
        CASE (15) ; CALL SHAPE_15_T (PT_L (1), PT_L (2), H_L)

!-->       QUADRILATERAL 2-D ELEMENTS
        CASE ( 4) ; CALL SHAPE_4_Q  (PT_L (1), PT_L (2), H_L)
        CASE ( 8) ; CALL SHAPE_8_Q  (PT_L (1), PT_L (2), H_L)
        CASE ( 9) ; CALL SHAPE_9_Q  (PT_L (1), PT_L (2), H_L)
        CASE (12) ; CALL SHAPE_12_Q (PT_L (1), PT_L (2), H_L)
        CASE (16) ; CALL SHAPE_16_Q (PT_L (1), PT_L (2), H_L)
        CASE (17) ; CALL SHAPE_17_Q (PT_L (1), PT_L (2), H_L)
!       CASE (25) ; CALL SHAPE_25_Q (PT_L (1), PT_L (2), H_L)
        CASE DEFAULT 
          PRINT *, 'WARNING: CHECK VALUE OF shape KEYWORD'
          N_WARN = N_WARN + 1
          CALL SHAPE_USER (PT_L, H_L, LT_N, LT_PARM, N_G_DOF)
      END SELECT ! number of element nodes
    CASE (3) ! 3-D space
      SELECT CASE ( LT_N )

!-->       HEXAHEDRA 3-D ELEMENTS
        CASE ( 8) ; CALL SHAPE_8_H  (PT_L (1), PT_L (2), PT_L (3), H_L)
        CASE (20) ; CALL SHAPE_20_H (PT_L (1), PT_L (2), PT_L (3), H_L)
!       CASE (27) ; CALL SHAPE_27_H (PT_L (1), PT_L (2), PT_L (3), H_L)
!       CASE (32) ; CALL SHAPE_32_H (PT_L (1), PT_L (2), PT_L (3), H_L)

!-->       TETRAHEDRA 3-D ELEMENTS (PYRAMIDS)
        CASE ( 4) ; CALL SHAPE_4_P  (PT_L (1), PT_L (2), PT_L (3), H_L)
!       CASE (10) ; CALL SHAPE_10_P (PT_L (1), PT_L (2), PT_L (3), H_L)
!       CASE (20) ; CALL SHAPE_20_P (PT_L (1), PT_L (2), PT_L (3), H_L)

!-->       WEDGE 3-D ELEMENTS
!       CASE ( 6) ; CALL SHAPE_6_W  (PT_L (1), PT_L (2), PT_L (3), H_L)
!       CASE (15) ; CALL SHAPE_15_W (PT_L (1), PT_L (2), PT_L (3), H_L)

!-->       USER SUPPLIED ELEMENT
        CASE DEFAULT 
          CALL SHAPE_USER (PT_L, H_L, LT_N, LT_PARM, N_G_DOF)
      END SELECT ! number of element nodes
    CASE DEFAULT
      PRINT *, LT_N, LT_PARM, N_G_DOF
      STOP 'INVALID SPACE SIZE IN SCALAR_SHAPES'
END SELECT ! number of spatial dimensions
END SUBROUTINE SCALAR_SHAPES
                         
SUBROUTINE ELEM_C_N_SHAPES (PT_L, LEN_1, H_L)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
! EVALUATE C_N ELEMENT INTERPOLATION FUNCTIONS IN 1D, 2D, 3D
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants
Use Elem_Type_Data
 IMPLICIT NONE
 REAL(DP), INTENT(IN)  :: PT_L (LT_PARM), LEN_1
 REAL(DP), INTENT(OUT) :: H_L  (LT_N) !b XXX must be bigger
                                                  
! H_L       = ELEMENT INTERPOLATION FUNCTIONS AT PT_L
! LEN_1   = PHYSICAL LENGTH IN PARAMETRIC DIRECTION 1
! LT_N    = NUMBER OF NODES PER ELEMENT
! N_G_DOF = NUMBER OF DEGREES OF FREEDOM PER NODE
! LT_PARM = NO OF SPATIAL DIMENSIONS
! PT_L      = LOCAL COORD OF A POINT

!  BRANCH ON SPACE, THEN NUMBER OF NODES 
  SELECT CASE ( LT_PARM )
    CASE (0) ! point
      PRINT *,'WARNING, NO POINT ELEMENT C_N SHAPES'     
      N_WARN = N_WARN + 1 ; H = 1.d0
    CASE (1) ! 1-D space

!-->     1-D ELEMENTS
      SELECT CASE ( LT_N )
        CASE (2) 
          SELECT CASE ( N_G_DOF )
            CASE (2) ; CALL SHAPE_C1_L (PT_L (1), LEN_1, H_L)
            CASE (3) ; CALL SHAPE_C2_L (PT_L (1), LEN_1, H_L)
            CASE (4) ; CALL SHAPE_C3_L (PT_L (1), LEN_1, H_L)
            CASE DEFAULT 
              CALL SHAPE_C_N_USER (PT_L, LEN_1, H_L, LT_N, &
                                  LT_PARM, N_G_DOF)
          END SELECT
        CASE DEFAULT 
          CALL SHAPE_C_N_USER (PT_L, LEN_1, H_L, LT_N, &
                              LT_PARM, N_G_DOF)
      END SELECT ! number of element nodes
    CASE (2) ! 2-D space
      SELECT CASE ( LT_N )

!-->       RECTANGULAR 2-D ELEMENTS
!       CASE (17) 
!         CALL SHAPE_16_R (PT_L (1), PT_L (2), LEN_1, LEN_2, H_L)
        CASE DEFAULT 
          CALL SHAPE_C_N_USER (PT_L, LEN_1, H_L, LT_N, &
                              LT_PARM, N_G_DOF)
      END SELECT ! number of element nodes
    CASE DEFAULT
      PRINT *, LT_N, LT_PARM, N_G_DOF
      STOP 'INVALID SPACE SIZE IN ELEM_C_N_SHAPES'
END SELECT ! number of spatial dimensions
END SUBROUTINE ELEM_C_N_SHAPES
                         
SUBROUTINE GEOM_DERIVS (PT_L, DLG_L)   
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
! C0 ELEMENT GEOMETRIC DERIVATIVE FUNCTIONS IN 1D, 2D, 3D
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants
Use Elem_Type_Data
 IMPLICIT NONE
 REAL(DP), INTENT(IN)  :: PT_L  (LT_PARM)
 REAL(DP), INTENT(OUT) :: DLG_L (LT_PARM, LT_GEOM)
                                                  
! DLG_L     = LOCAL DERIVATIVES OF GEOMETRIC INTERPOLATIONS
! LT_GEOM = NUMBER OF GEOMETRIC NODES PER ELEMENT
! LT_PARM = NO OF SPATIAL DIMENSIONS
! PT_L      = LOCAL COORD OF A POINT
!
!  BRANCH ON SPACE, THEN NUMBER OF NODES 
  SELECT CASE ( LT_PARM )
    CASE (0) ! point
      PRINT *,'WARNING, NO POINT GEOMETRIC DERIVATIVES'
      N_WARN = N_WARN + 1 ; 
    CASE (1) ! 1-D space

!-->     1-D ELEMENTS
      SELECT CASE ( LT_GEOM )
        CASE (2) ; CALL DERIV_2_L (PT_L (1), DLG_L)
        CASE (3) ; CALL DERIV_3_L (PT_L (1), DLG_L)
        CASE (4) ; CALL DERIV_4_L (PT_L (1), DLG_L)
        CASE DEFAULT 
          CALL DERIV_USER (PT_L, DLG_L, LT_GEOM, LT_PARM, 1)      
      END SELECT ! number of element nodes
    CASE (2) ! 2-D space
      SELECT CASE ( LT_GEOM )

!-->       LINE ELEMENT IN 2_D
        CASE ( 2) ; CALL DERIV_2_L  (PT_L (1), DLG_L)

!-->       TRIANGULAR 2-D ELEMENTS
        CASE ( 3) ; CALL DERIV_3_T  (PT_L (1), PT_L (2), DLG_L)
!       CASE ( 4) ; CALL DERIV_4_T  (PT_L (1), PT_L (2), DLG_L)
        CASE ( 6) ; CALL DERIV_6_T  (PT_L (1), PT_L (2), DLG_L)
!       CASE ( 7) ; CALL DERIV_7_T  (PT_L (1), PT_L (2), DLG_L)
        CASE (10) ; CALL DERIV_10_T (PT_L (1), PT_L (2), DLG_L)
        CASE (15) ; CALL DERIV_15_T (PT_L (1), PT_L (2), DLG_L)

!-->       QUADRILATERAL 2-D ELEMENTS
        CASE ( 4) ; CALL DERIV_4_Q  (PT_L (1), PT_L (2), DLG_L)
        CASE ( 8) ; CALL DERIV_8_Q  (PT_L (1), PT_L (2), DLG_L)
        CASE ( 9) ; CALL DERIV_9_Q  (PT_L (1), PT_L (2), DLG_L)
        CASE (12) ; CALL DERIV_12_Q (PT_L (1), PT_L (2), DLG_L)
        CASE (16) ; CALL DERIV_16_Q (PT_L (1), PT_L (2), DLG_L)
        CASE (17) ; CALL DERIV_17_Q (PT_L (1), PT_L (2), DLG_L)
!       CASE (25) ; CALL DERIV_25_Q (PT_L (1), PT_L (2), DLG_L)
        CASE DEFAULT 
          CALL DERIV_USER (PT_L, DLG_L, LT_GEOM, LT_PARM, 1)      
      END SELECT ! number of element nodes
    CASE (3) ! 3-D space
      SELECT CASE ( LT_GEOM )

!-->       LINE ELEMENT IN 3_D
        CASE ( 2) ; CALL DERIV_2_L  (PT_L (1), DLG_L)

!-->       TRIANGULAR 2-D ELEMENTS IN 3-D
        CASE ( 3) ; CALL DERIV_3_T  (PT_L (1), PT_L (2), DLG_L)
        CASE ( 6) ; CALL DERIV_6_T  (PT_L (1), PT_L (2), DLG_L)

!-->       QUADRILATERAL 2-D ELEMENTS IN 3-D, OR TET
        CASE ( 4) 
          IF ( LT_SHAP == 3 ) CALL DERIV_4_Q (PT_L (1), PT_L (2), DLG_L)
          IF ( LT_SHAP == 5 ) CALL DERIV_4_P (PT_L (1), PT_L (2), &
                                              PT_L (3), DLG_L)
        CASE ( 9) ; CALL DERIV_9_Q  (PT_L (1), PT_L (2), DLG_L)
        CASE (12) ; CALL DERIV_12_Q (PT_L (1), PT_L (2), DLG_L)
        CASE (16) ; CALL DERIV_16_Q (PT_L (1), PT_L (2), DLG_L)

!-->       HEXAHEDRA 3-D ELEMENTS
        CASE ( 8) ; CALL DERIV_8_H  (PT_L (1), PT_L (2), PT_L (3), DLG_L)
        CASE (20) ; CALL DERIV_20_H (PT_L (1), PT_L (2), PT_L (3), DLG_L)
!       CASE (27) ; CALL DERIV_27_H (PT_L (1), PT_L (2), PT_L (3), DLG_L)
!       CASE (32) ; CALL DERIV_32_H (PT_L (1), PT_L (2), PT_L (3), DLG_L)

!-->       TETRAHEDRA 3-D ELEMENTS (PYRAMIDS)
!       CASE (10) ; CALL DERIV_10_P (PT_L (1), PT_L (2), PT_L (3), DLG_L)
!       CASE (20) ; CALL DERIV_20_P (PT_L (1), PT_L (2), PT_L (3), DLG_L)

!-->       WEDGE 3-D ELEMENTS
!       CASE ( 6) ; CALL DERIV_6_W  (PT_L (1), PT_L (2), PT_L (3), DLG_L)
!       CASE (15) ; CALL DERIV_15_W (PT_L (1), PT_L (2), PT_L (3), DLG_L)

!-->       USER SUPPLIED ELEMENT
        CASE DEFAULT 
          CALL DERIV_USER (PT_L, DLG_L, LT_GEOM, LT_PARM, 1)      
      END SELECT ! number of element nodes
    CASE DEFAULT
      PRINT *, LT_GEOM, LT_PARM
      STOP 'INVALID SPACE SIZE IN GEOM_DERIVS'
END SELECT ! number of spatial dimensions
END SUBROUTINE GEOM_DERIVS
                         
SUBROUTINE GEOM_SHAPES (PT_L, G_L)   
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
! C0 ELEMENT GEOMETRIC INTERPOLATION FUNCTIONS IN 1D, 2D, 3D
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants
Use Elem_Type_Data
 IMPLICIT NONE
 REAL(DP), INTENT(IN)  :: PT_L (LT_PARM)
 REAL(DP), INTENT(OUT) :: G_L  (LT_GEOM)
                                                  
! G_L       = ELEMENT INTERPOLATION FUNCTIONS AT PT_L
! LT_GEOM = NUMBER OF GEOMETRIC NODES PER ELEMENT
! LT_PARM = NO OF SPATIAL DIMENSIONS
! PT_L      = LOCAL COORD OF A POINT

!  BRANCH ON SPACE, THEN NUMBER OF NODES 
  SELECT CASE ( LT_PARM )
    CASE (0) ! point
      PRINT *,'WARNING, NO POINT GEOMETRIC SHAPES'     
      N_WARN = N_WARN + 1 ; G_L = 1.d0
    CASE (1) ! 1-D space

!-->     1-D ELEMENTS
      SELECT CASE ( LT_GEOM )
        CASE (2) ; CALL SHAPE_2_L (PT_L (1), G_L)
        CASE (3) ; CALL SHAPE_3_L (PT_L (1), G_L)
        CASE (4) ; CALL SHAPE_4_L (PT_L (1), G_L)
        CASE DEFAULT 
          CALL SHAPE_USER (PT_L, G_L, LT_GEOM, LT_PARM, 1)      
      END SELECT ! number of element nodes
    CASE (2) ! 2-D space
      SELECT CASE ( LT_GEOM )

!-->       LINE ELEMENT IN 2-D
        CASE (2) ; CALL SHAPE_2_L (PT_L (1), G_L)

!-->       TRIANGULAR 2-D ELEMENTS
        CASE ( 3) ; CALL SHAPE_3_T  (PT_L (1), PT_L (2), G_L)
!       CASE ( 4) ; CALL SHAPE_4_T  (PT_L (1), PT_L (2), G_L)
        CASE ( 6) ; CALL SHAPE_6_T  (PT_L (1), PT_L (2), G_L)
!       CASE ( 7) ; CALL SHAPE_7_T  (PT_L (1), PT_L (2), G_L)
        CASE (10) ; CALL SHAPE_10_T (PT_L (1), PT_L (2), G_L)
        CASE (15) ; CALL SHAPE_15_T (PT_L (1), PT_L (2), G_L)

!-->       QUADRILATERAL 2-D ELEMENTS
        CASE ( 4) ; CALL SHAPE_4_Q  (PT_L (1), PT_L (2), G_L)
        CASE ( 8) ; CALL SHAPE_8_Q  (PT_L (1), PT_L (2), G_L)
        CASE ( 9) ; CALL SHAPE_9_Q  (PT_L (1), PT_L (2), G_L)
        CASE (12) ; CALL SHAPE_12_Q (PT_L (1), PT_L (2), G_L)
        CASE (16) ; CALL SHAPE_16_Q (PT_L (1), PT_L (2), G_L)
        CASE (17) ; CALL SHAPE_17_Q (PT_L (1), PT_L (2), G_L)
!       CASE (25) ; CALL SHAPE_25_Q (PT_L (1), PT_L (2), G_L)
        CASE DEFAULT 
          CALL SHAPE_USER (PT_L, G_L, LT_GEOM, LT_PARM, 1)      
      END SELECT ! number of element nodes
    CASE (3) ! 3-D space
      SELECT CASE ( LT_GEOM )

!-->       LINE ELEMENT IN 3-D
        CASE (2) ; CALL SHAPE_2_L (PT_L (1), G_L)

!-->       TRIANGULAR 2-D ELEMENTS IN 3-D
        CASE ( 3) ; CALL SHAPE_3_T  (PT_L (1), PT_L (2), G_L)
        CASE ( 6) ; CALL SHAPE_6_T  (PT_L (1), PT_L (2), G_L)

!-->       QUADRILATERAL 2-D ELEMENTS IN 3-D, OR TET
        CASE ( 4) 
          IF ( LT_SHAP == 3 ) CALL SHAPE_4_Q (PT_L (1), PT_L (2), G_L)
          IF ( LT_SHAP == 5 ) CALL SHAPE_4_P (PT_L (1), PT_L (2),    &
                                              PT_L (3), G_L)
        CASE ( 9) ; CALL SHAPE_9_Q  (PT_L (1), PT_L (2), G_L)
        CASE (12) ; CALL SHAPE_12_Q (PT_L (1), PT_L (2), G_L)
        CASE (16) ; CALL SHAPE_16_Q (PT_L (1), PT_L (2), G_L)

!-->       HEXAHEDRA 3-D ELEMENTS
        CASE ( 8) ; CALL SHAPE_8_H  (PT_L (1), PT_L (2), PT_L (3), G_L)
        CASE (20) ; CALL SHAPE_20_H (PT_L (1), PT_L (2), PT_L (3), G_L)
!       CASE (27) ; CALL SHAPE_27_H (PT_L (1), PT_L (2), PT_L (3), G_L)
!       CASE (32) ; CALL SHAPE_32_H (PT_L (1), PT_L (2), PT_L (3), G_L)

!-->       TETRAHEDRA 3-D ELEMENTS (PYRAMIDS)
!b      CASE ( 4) ; CALL SHAPE_4_P  (PT_L (1), PT_L (2), PT_L (3), G_L)
!       CASE (10) ; CALL SHAPE_10_P (PT_L (1), PT_L (2), PT_L (3), G_L)
!       CASE (20) ; CALL SHAPE_20_P (PT_L (1), PT_L (2), PT_L (3), G_L)

!-->       WEDGE 3-D ELEMENTS
!       CASE ( 6) ; CALL SHAPE_6_W  (PT_L (1), PT_L (2), PT_L (3), G_L)
!       CASE (15) ; CALL SHAPE_15_W (PT_L (1), PT_L (2), PT_L (3), G_L)

!-->       USER SUPPLIED ELEMENT
        CASE DEFAULT 
          CALL SHAPE_USER (PT_L, G_L, LT_GEOM, LT_PARM, 1)      
      END SELECT ! number of element nodes
    CASE DEFAULT
      PRINT *, LT_GEOM, LT_PARM
      STOP 'INVALID SPACE SIZE IN GEOM_SHAPES'
END SELECT ! number of spatial dimensions
END SUBROUTINE GEOM_SHAPES
                         
SUBROUTINE GEN_ELEM_SHAPE (PT_L, H, NOD_PER_EL, N_SPACE, N_G_DOF)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
! EVALUATE C0 ELEMENT INTERPOLATION FUNCTIONS IN 1D, 2D, 3D
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
Use Precision_Module
 IMPLICIT NONE
 INTEGER,  INTENT(IN)  :: NOD_PER_EL, N_SPACE, N_G_DOF
 REAL(DP), INTENT(IN)  :: PT_L (N_SPACE)
 REAL(DP), INTENT(OUT) :: H    (NOD_PER_EL * N_G_DOF)
                                                  
! H          = ELEMENT INTERPOLATION FUNCTIONS AT PT_L
! NOD_PER_EL = NUMBER OF NODES PER ELEMENT
! N_G_DOF    = NUMBER OF DEGREES OF FREEDOM PER NODE
! N_SPACE    = NO OF SPATIAL DIMENSIONS
! PT_L         = LOCAL COORD OF A POINT

!  BRANCH ON SPACE, THEN NUMBER OF NODES (and then N_G_DOF)
  SELECT CASE ( N_SPACE )
    CASE (1) ! 1-D space

!-->     1-D ELEMENTS
      SELECT CASE ( NOD_PER_EL )
        CASE (2) ; CALL SHAPE_2_L (PT_L (1), H)
        CASE (3) ; CALL SHAPE_3_L (PT_L (1), H)
        CASE (4) ; CALL SHAPE_4_L (PT_L (1), H)
        CASE DEFAULT 
          CALL SHAPE_USER (PT_L, H, NOD_PER_EL, N_SPACE, N_G_DOF)
      END SELECT ! number of element nodes
    CASE (2) ! 2-D space
      SELECT CASE ( NOD_PER_EL )

!-->       LINE ELEMENTS IN 2-D
        CASE ( 2) ; CALL SHAPE_2_L  (PT_L (1), H)

!-->       TRIANGULAR 2-D ELEMENTS
        CASE ( 3) ; CALL SHAPE_3_T  (PT_L (1), PT_L (2), H)
!       CASE ( 4) ; CALL SHAPE_4_T  (PT_L (1), PT_L (2), H)
        CASE ( 6) ; CALL SHAPE_6_T  (PT_L (1), PT_L (2), H)
!       CASE ( 7) ; CALL SHAPE_7_T  (PT_L (1), PT_L (2), H)
        CASE (10) ; CALL SHAPE_10_T (PT_L (1), PT_L (2), H)
        CASE (15) ; CALL SHAPE_15_T (PT_L (1), PT_L (2), H)

!-->       QUADRILATERAL 2-D ELEMENTS
        CASE ( 4) ; CALL SHAPE_4_Q  (PT_L (1), PT_L (2), H)
        CASE ( 8) ; CALL SHAPE_8_Q  (PT_L (1), PT_L (2), H)
        CASE ( 9) ; CALL SHAPE_9_Q  (PT_L (1), PT_L (2), H)
        CASE (12) ; CALL SHAPE_12_Q (PT_L (1), PT_L (2), H)
        CASE (16) ; CALL SHAPE_16_Q (PT_L (1), PT_L (2), H)
        CASE (17) ; CALL SHAPE_17_Q (PT_L (1), PT_L (2), H)
!       CASE (25) ; CALL SHAPE_25_Q (PT_L (1), PT_L (2), H)
        CASE DEFAULT 
          CALL SHAPE_USER (PT_L, H, NOD_PER_EL, N_SPACE, N_G_DOF)
      END SELECT ! number of element nodes
    CASE (3) ! 3-D space
      SELECT CASE ( NOD_PER_EL )

!-->       LINE ELEMENTS IN 3-D
        CASE ( 2) ; CALL SHAPE_2_L  (PT_L (1), H)

!-->       HEXAHEDRA 3-D ELEMENTS
        CASE ( 8) ; CALL SHAPE_8_H  (PT_L (1), PT_L (2), PT_L (3), H)
        CASE (20) ; CALL SHAPE_20_H (PT_L (1), PT_L (2), PT_L (3), H)
!       CASE (27) ; CALL SHAPE_27_H (PT_L (1), PT_L (2), PT_L (3), H)
!       CASE (32) ; CALL SHAPE_32_H (PT_L (1), PT_L (2), PT_L (3), H)

!-->       TETRAHEDRA 3-D ELEMENTS (PYRAMIDS)
        CASE ( 4) ; CALL SHAPE_4_P  (PT_L (1), PT_L (2), PT_L (3), H)
!       CASE (10) ; CALL SHAPE_10_P (PT_L (1), PT_L (2), PT_L (3), H)
!       CASE (20) ; CALL SHAPE_20_P (PT_L (1), PT_L (2), PT_L (3), H)

!-->       WEDGE 3-D ELEMENTS
!       CASE ( 6) ; CALL SHAPE_6_W  (PT_L (1), PT_L (2), PT_L (3), H)
!       CASE (15) ; CALL SHAPE_15_W (PT_L (1), PT_L (2), PT_L (3), H)

!-->       USER SUPPLIED ELEMENT
        CASE DEFAULT 
          CALL SHAPE_USER (PT_L, H, NOD_PER_EL, N_SPACE, N_G_DOF)
      END SELECT ! number of element nodes
    CASE DEFAULT
      STOP 'INVALID SPACE SIZE IN GEN_ELEM_SHAPE'
END SELECT ! number of spatial dimensions
END SUBROUTINE GEN_ELEM_SHAPE

SUBROUTINE DERIV_10_T (R, S, DH)
! **********************************************************
! LOCAL DERIVATIVES OF SHAPE FUNCTIONS FOR 10-NODED TRIANGLE
! **********************************************************
Use Precision_Module
 IMPLICIT NONE
 REAL(DP), INTENT(IN)  :: R, S
 REAL(DP), INTENT(OUT) :: DH (2, 10)
 REAL(DP) :: T0, T1, T2, T3, T4, T5, T6, T7, T8, T9, T10, Z

!                                   S
! THE CUBIC TRIANGLE                :
!                                   3
! R,S = LOCAL COORDS OF RT          | \
! H = ELEM SHARE FUNCTIONS          |   \
!                                   6     8
! ELEM_NODES = ELEM INCIDENCES LIST |      \
! ELEMENT SKETCH TO RIGHT           9   10   5
!                                   |          \
! 1@(0,0)  2@(1,0) 2@(0,1)          1---4---7---2  ... R

  Z   = 1.d0 - R - S
  T0  = 4.5d0 * (6.d0 * Z - 1.d0) 
  T1  = 4.5d0 * S * (3.d0 * S - 1.d0)
  T2  = 4.5d0 * Z * (3.d0 * Z - 1.d0)
  T3  = 4.5d0 * R * (3.d0 * R - 1.d0)
  T4  = 4.5d0 * (6.d0 * S - 1.d0)
  T5  = 4.5d0 * (6.d0 * R - 1.d0)
  T6  = T0 * S
  T7  = T0 * R
  T8  = -0.5d0 * (27.d0 * Z * Z - 18.d0 * Z + 2.d0)
  T9  =  0.5d0 * (27.d0 * S * S - 18.d0 * S + 2.d0)
  T10 =  0.5d0 * (27.d0 * R * R - 18.d0 * R + 2.d0)

  DH (1, 1 ) = T8
  DH (1, 2 ) = T10
  DH (1, 3 ) = 0.d0
  DH (1, 4 ) = T2 - T7
  DH (1, 5 ) = S * T5
  DH (1, 6 ) = -T1
  DH (1, 7 ) = Z * T5 - T3
  DH (1, 8 ) = T1
  DH (1, 9 ) = -T6
  DH (1, 10) = 27.d0 * S * (Z - R)

  DH (2, 1 ) = T8
  DH (2, 2 ) = 0.d0
  DH (2, 3 ) = T9
  DH (2, 4 ) = -T7
  DH (2, 5 ) = T3
  DH (2, 6 ) = Z * T4 - T1
  DH (2, 7 ) = -T3
  DH (2, 8 ) = R * T4
  DH (2, 9 ) = T2 - T6
  DH (2, 10) = 27.d0 * R * (Z - S)
END SUBROUTINE DERIV_10_T 

SUBROUTINE SHAPE_10_T (R, S, H)
! ******************************************************
! SHAPE FUNCTIONS FOR 10-NODED TRIANGLE
! ******************************************************
Use Precision_Module
 IMPLICIT NONE
 REAL(DP), INTENT(IN)  :: R, S
 REAL(DP), INTENT(OUT) :: H (10)
 REAL(DP) :: T1, T2, T3, Z

!                                   S
! THE CUBIC TRIANGLE                :
!                                   3
! R,S = LOCAL COORDS OF PT          | \
! H = ELEM SHAPE FUNCTIONS          |   \
!                                   6     8
! ELEM_NODES = ELEM INCIDENCES LIST |      \
! ELEMENT SKETCH TO RIGHT           9   10   5
!                                   |          \
! 1@(0,0)  2@(1,0) 2@(0,1)          1---4---7---2  ... R

  Z  = 1.d0 - R - S
  T1 = S * (3.d0 * S - 1.d0)
  T2 = Z * (3.d0 * Z - 1.d0)
  T3 = R * (3.d0 * R - 1.d0)
  H (1 ) = 0.5d0 * T2 * (3.d0 * Z - 2.d0)
  H (2 ) = 0.5d0 * T3 * (3.d0 * R - 2.d0)
  H (3 ) = 0.5d0 * T1 * (3.d0 * S - 2.d0)
  H (4 ) = 4.5d0 * R * T2
  H (5 ) = 4.5d0 * S * T3
  H (6 ) = 4.5d0 * Z * T1
  H (7 ) = 4.5d0 * Z * T3
  H (8 ) = 4.5d0 * R * T1
  H (9 ) = 4.5d0 * S * T2
  H (10) = 27.d0 * S * Z * R
END SUBROUTINE SHAPE_10_T
                         
SUBROUTINE DERIV_12_Q (R, S, DH)
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *
! DERIVATIVES OF A 12 NODE QUAD, SEE SHAPE_12_Q
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *
Use Precision_Module
 IMPLICIT NONE
 REAL(DP), INTENT(IN)  :: R, S
 REAL(DP), INTENT(OUT) :: DH (2, 12)
 REAL(DP), PARAMETER   :: F1 = 1.d0/32, F2 = 9.d0/32
 REAL(DP) :: R_P, R_M, S_P, S_M, R2, S2

! R,S = LOCAL COORDS OF PT     4--11---7---3
! H = ELEM SHAPE FUNCTIONS     I           I
!                              8     S    10
! ELEMENT SKETCH TO RIGHT      I     .R    I
! LOCAL COORD OF NODES:       12           6
! 1-(-1,-1)  3-(+1,+1)         I           I
! SIDES OF ORDER 1, 2, OR 3    1---5---9---2

  R_P = 1 + R ; R_M = 1 - R ; S_P = 1 + S ; S_M = 1 - S
  R2 = R*R ; S2 = S*S
  DH (1,  1) =  F1*S_M*( 10 - 27*R2 + 18*R - 9*S2)
  DH (2,  1) =  F1*R_M*( 10 - 27*S2 + 18*S - 9*R2)
  DH (1,  2) =  F1*S_M*(-10 + 27*R2 + 18*R + 9*S2)
  DH (2,  2) =  F1*R_P*( 10 - 27*S2 + 18*S - 9*R2)
  DH (1,  3) =  F1*S_P*(-10 + 27*R2 + 18*R + 9*S2)
  DH (2,  3) =  F1*R_P*(-10 + 27*S2 + 18*S + 9*R2)
  DH (1,  4) =  F1*S_P*( 10 - 27*R2 + 18*R - 9*S2)
  DH (2,  4) =  F1*R_M*(-10 + 27*S2 + 18*S + 9*R2)
  DH (1,  5) =  F2*S_M*( 9*R2 - 2*R - 3)
  DH (2,  5) = -F2*(1 - R2)*(1 - 3*R)
  DH (1,  6) =  F2*(1 - S2)*(1 - 3*S)
  DH (2,  6) =  F2*R_P*( 9*S2 - 2*S - 3)
  DH (1,  7) =  F2*S_P*(-9*R2 - 2*R + 3)
  DH (2,  7) =  F2*(1 - R2)*(1 + 3*R)
  DH (1,  8) = -F2*(1 - S2)*(1 + 3*S)
  DH (2,  8) =  F2*R_M*(-9*S2 - 2*S + 3)
  DH (1,  9) =  F2*S_M*(-9*R2 - 2*R + 3)
  DH (2,  9) = -F2*(1 - R2)*(1 + 3*R)
  DH (1, 10) =  F2*(1 - S2)*(1 + 3*S)
  DH (2, 10) =  F2*R_P*(-9*S2 - 2*S + 3)
  DH (1, 11) =  F2*S_P*( 9*R2 - 2*R - 3)
  DH (2, 11) =  F2*(1 - R2)*(1 - 3*R)
  DH (1, 12) = -F2*(1 - S2)*(1 - 3*S)
  DH (2, 12) =  F2*R_M*( 9*S2 - 2*S - 3)
END SUBROUTINE DERIV_12_Q
                         
SUBROUTINE DERIV_15_T (R, S, DH)
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-*
! LOCAL DERIVATIVES FOR 15-NODED TRIANGLE
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-*
Use Precision_Module
 IMPLICIT NONE
 REAL(DP), INTENT(IN)  :: R, S
 REAL(DP), PARAMETER   :: PT1 = 42.666666666666667d0, &
                          PT2 = 10.666666666666667d0
 REAL(DP), INTENT(OUT) :: DH (2, 15)
 REAL(DP) :: CC, T1, T2, T3, T4, T5, T6, T7, T8, T9
  CC = 1.d0 - R - S
  T1 = R - 0.25d0  ; T2 = R - 0.5d0  ; T3 = R - 0.75d0
  T4 = S - 0.25d0  ; T5 = S - 0.5d0  ; T6 = S - 0.75d0
  T7 = CC - 0.25d0 ; T8 = CC - 0.5d0 ; T9 = CC - 0.75d0
  DH (1, 2) = PT2 * (T2 * T3 * (T1 + R) + R * T1 * (T3 + T2) )
  DH (1, 5) = PT1 * S * (T2 * (T1 + R) + R * T1)
  DH (1, 8) = 64.d0 * S * T4 * (T1 + R)
  DH (1, 11) = PT1 * S * T4 * T5
  DH (1, 3) = 0.d0
  DH (1, 6) = - PT1 * S * T4 * T5
  DH (1, 9) = - 64.d0 * S * T4 * (T7 + CC)
  DH (1, 12) = - PT1 * S * (T8 * (T7 + CC) + CC * T7)
  DH (1, 1) = - PT2 * (T8 * T9 * (T7 + CC) + CC * T7 * (T8 + T9) )
  DH (1, 4) = PT1 * (CC * T7 * T8 - R * (T8 * (T7 + CC) + CC * T7) )
  DH (1, 7) = 64.d0 * (CC * T7 * (T1 + R) - R * T1 * (T7 + CC) )
  DH (1, 10) = PT1 * (CC * (T2 * (T1 + R) + R * T1) - R * T1 * T2)
  DH (1, 14) = 128.d0 * S * (CC * (T1 + R) - R * T1)
  DH (1, 15) = 128.d0 * S * T4 * (CC - R)
  DH (1, 13) = 128.d0 * S * (CC * T7 - R * (T7 + CC) )
  DH (2, 2) = 0.d0
  DH (2, 5) = PT1 * R * T1 * T2
  DH (2, 8) = 64.d0 * R * T1 * (T4 + S)
  DH (2, 11) = PT1 * R * (T5 * (T4 + S) + S * T4)
  DH (2, 3) = PT2 * (T5 * T6 * (T4 + S) + S * T4 * (T6 + T5) )
  DH (2, 6) = PT1 * ((CC * (T5 * (T4 + S) + S * T4)) - S * T4 * T5)
  DH (2, 9) = 64.d0 * (CC * T7 * (T4 + S) - S * T4 * (T7 + CC) )
  DH (2, 12) = PT1 * (CC * T7 * T8 - S * (T8 * (T7 + CC) + CC * T7))
  DH (2, 1) = - PT2 * (T8 * T9 * (T7 + CC) + CC * T7 * (T8 + T9) )
  DH (2, 4) = - PT1 * R * (T8 * (T7 + CC) + CC * T7)
  DH (2, 7) = - 64.d0 * R * T1 * (T7 + CC)
  DH (2, 10) = - PT1 * R * T1 * T2
  DH (2, 14) = 128.d0 * R * T1 * (CC - S)
  DH (2, 15) = 128.d0 * R * (CC * (T4 + S) - S * T4)
  DH (2, 13) = 128.d0 * R * (CC * T7 - S * (CC + T7) )
END SUBROUTINE DERIV_15_T

SUBROUTINE DERIV_16_Q (R, S, DH)
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *
! DERIVATIVES OF A LAGRANGIAN 16 NODE QUAD
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *
!   SEE SHAPE_16_Q
Use Precision_Module
 IMPLICIT NONE
 REAL(DP), INTENT(IN)  :: R, S
 REAL(DP), INTENT(OUT) :: DH (2, 16)
 REAL(DP)              :: HR (4), HS (4), DHR (4), DHS (4)

! R,S = LOCAL COORDS OF PT     4--11---7---3
! H = ELEM SHAPE FUNCTIONS     I           I
!                              8  16 S 15 10
! ELEMENT SKETCH TO RIGHT      I     .R    I
! LOCAL COORD OF NODES:       12  13   14  6
! 1-(-1,-1)  3-(+1,+1)         I           I
!                              1---5---9---2

! Evaluate the 1D interpolations
  CALL SHAPE_4_L (R, HR)
  CALL SHAPE_4_L (S, HS)
  CALL DERIV_4_L (R, DHR)
  CALL DERIV_4_L (S, DHS)

  DH (1,  1) =  DHR (1)*HS (1)
  DH (2,  1) =  HR (1)*DHS (1)
  DH (1,  2) =  DHR (4)*HS (1)
  DH (2,  2) =  HR (4)*DHS (1)
  DH (1,  3) =  DHR (4)*HS (4)
  DH (2,  3) =  HR (4)*DHS (4)
  DH (1,  4) =  DHR (1)*HS (4)
  DH (2,  4) =  HR (1)*DHS (4)
  DH (1,  5) =  DHR (2)*HS (1)
  DH (2,  5) =  HR (2)*DHS (1)
  DH (1,  6) =  DHR (4)*HS (2)
  DH (2,  6) =  HR (4)*DHS (2)
  DH (1,  7) =  DHR (3)*HS (4)
  DH (2,  7) =  HR (3)*DHS (4)
  DH (1,  8) =  DHR (1)*HS (3)
  DH (2,  8) =  HR (1)*DHS (3)
  DH (1,  9) =  DHR (3)*HS (1)
  DH (2,  9) =  HR (3)*DHS (1)
  DH (1, 10) =  DHR (4)*HS (3)
  DH (2, 10) =  HR (4)*DHS (3)
  DH (1, 11) =  DHR (2)*HS (4)
  DH (2, 11) =  HR (2)*DHS (4)
  DH (1, 12) =  DHR (1)*HS (2)
  DH (2, 12) =  HR (1)*DHS (2)
  DH (1, 13) =  DHR (2)*HS (2)
  DH (2, 13) =  HR (2)*DHS (2)
  DH (1, 14) =  DHR (3)*HS (2)
  DH (2, 14) =  HR (3)*DHS (2)
  DH (1, 15) =  DHR (3)*HS (3)
  DH (2, 15) =  HR (3)*DHS (3)
  DH (1, 16) =  DHR (2)*HS (3)
  DH (2, 16) =  HR (2)*DHS (3)
END SUBROUTINE DERIV_16_Q
                         
SUBROUTINE DERIV_16_QS (R, S, DH)
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *
! LOCAL DERIVATIVES FOR A SERENDIPITY 16 NODE QUAD
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *
! SEE SHAPE_16_QS FOR NODE LOCATIONS
Use Precision_Module
 IMPLICIT NONE
 REAL(DP), INTENT(IN)  :: R, S
 REAL(DP), PARAMETER   :: PT667  = 0.66666666666667d0, &
                          PT1D12 = 0.083333333333333d0, HALF = 0.5d0
 REAL(DP), INTENT(OUT) :: DH (2, 16)
 REAL(DP) :: RR, SS, RRR, SSS, RP, RM, S_P, SM, R1
  RR = R * R       ; SS = S * S
  RRR = R * R * R  ; SSS = S * S * S
  RP = 1.d0 + R      ; RM = 1.d0 - R
  S_P = 1.d0 + S      ; SM = 1.d0 - S
  R1 = R - 1.d0
  DH (1, 1) = PT1D12*SM*(16.d0*RRR - 12.d0*RR - 2.d0*R + 4.d0*SSS - S + 4.d0)
  DH (1, 2) = PT1D12*SM*(16.d0*RRR + 12.d0*RR - 2.d0*R - 4.d0*SSS + S - 4.d0)
  DH (1, 3) = PT1D12*S_P*(16.d0*RRR + 12.d0*RR - 2.d0*R + 4.d0*SSS - S - 4.d0)
  DH (1, 4) = PT1D12*S_P*(16.d0*RRR - 12.d0*RR - 2.d0*R - 4.d0*SSS + S + 4.d0)
  DH (1, 5) = PT667 * (4.d0 * R - 1.d0 + 3.d0 * RR - 8.d0 * RRR) * SM
  DH (1, 6) = PT667 * S * SM * S_P * ( -1.d0 + 2.d0 * S)
  DH (1, 7) = PT667 * (1.d0 + 4.d0 * R - 3.d0 * RR - 8.d0 * RRR) * S_P
  DH (1, 8) = -PT667 * S * SM * S_P * (1.d0 + 2.d0 * S)
  DH (1, 9) = R * ( -5.d0 + 8.d0 * RR) * SM
  DH (1, 10) = HALF * SM * S_P * (1.d0 - 2.d0 * S) * (1.d0 + 2.d0 * S)
  DH (1, 11) = R * ( -5.d0 + 8.d0 * RR) * S_P
  DH (1, 12) = HALF * SM * S_P * ( -1.d0 + 2.d0 * S) * (1.d0 + 2.d0 * S)
  DH (1, 13) = PT667 * (1.d0 + 4.d0 * R - 3.d0 * RR - 8.d0 * RRR) * SM
  DH (1, 14) = PT667 * S * SM * S_P * (1.d0 + 2.d0 * S)
  DH (1, 15) = PT667 * (4.d0 * R - 1.d0 + 3.d0 * RR - 8.d0 * RRR) * S_P
  DH (1, 16) = PT667 * S * SM * S_P * (1.d0 - 2.d0 * S)
  DH (2, 1) = PT1D12 * RM * (16.d0 * SSS - 12.d0 * SS - 2.d0 * S &
                             + 4.d0 * RRR - R + 4.d0)
  DH (2, 2) = PT1D12 * RP * (16.d0 * SSS - 12.d0 * SS - 2.d0 * S &
                             - 4.d0 * RRR + R + 4.d0)
  DH (2, 3) = PT1D12 * RP * (16.d0 * SSS + 12.d0 * SS - 2.d0 * S &
                             + 4.d0 * RRR - R - 4.d0)
  DH (2, 4) = PT1D12 * R1 * ( -16.d0 * SSS - 12.d0 * SS &
                             + 2.d0 * S + 4.d0 * RRR - R + 4.d0)
  DH (2, 5) = PT667 * R * R1 * RP * (2.d0 * R - 1.d0)
  DH (2, 6) = PT667 * ( -1.d0 + 4.d0 * S + 3.d0 * SS - 8.d0 * SSS) * RP
  DH (2, 7) = PT667 * R * RM * RP * (1.d0 + 2.d0 * R)
  DH (2, 8) = PT667 * ( -1.d0 - 4.d0 * S + 3.d0 * SS + 8.d0 * SSS) * R1
  DH (2, 9) = HALF * RM * RP * (2.d0 * R - 1.d0) * (1.d0 + 2.d0 * R)
  DH (2, 10) = S * ( -5.d0 + 8.d0 * SS) * RP
  DH (2, 11) = HALF * R1 * RP * (2.d0 * R - 1.d0) * (1.d0 + 2.d0 * R)
  DH (2, 12) = S * (5.d0 - 8.d0 * SS) * R1
  DH (2, 13) = PT667 * R * R1 * RP * (1.d0 + 2.d0 * R)
  DH (2, 14) = PT667 * (1.d0 + 4.d0 * S - 3.d0 * SS - 8.d0 * SSS) * RP
  DH (2, 15) = PT667 * R * RM * RP * (2.d0 * R - 1.d0)
  DH (2, 16) = PT667 * (1.d0 - 4.d0 * S - 3.d0 * SS + 8.d0 * SSS) * R1
END SUBROUTINE DERIV_16_QS
                          
SUBROUTINE DERIV_16_R (R, S, A, B, DH)
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *
! WARNING, this code is not fully checked. No known bugs.
! FIRST DERIVATIVES OF A C1 RECTANGULAR ELEMENT IN
! UNIT COORDINATES USING TENSOR PRODUCTS OF 1D BASIS
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *
Use Precision_Module
 IMPLICIT NONE
 REAL(DP), INTENT(IN)  :: R, S, A, B
 REAL(DP), INTENT(OUT) :: DH (2, 16)
 REAL(DP)              :: HR (4), DHR (4), HS (4), DHS (4)

! DOF ARE W W,X W,Y W,XY AT EACH NODE (N_G_DOF=4)
! X // R, Y // S.              S
! A = PHYSICAL LENGTH IN X     4 -------- 3
! B = PHYSICAL LENGTH IN Y     I          I
! R,S = LOCAL UNIT COORDS      I          I
! 1@(0,0), 3@(1,1)             1 -------- 2 ->R

! Evaluate the 1D interpolations
  CALL SHAPE_C1_L (R, A, HR)
  CALL SHAPE_C1_L (S, B, HS)
  CALL DERIV_C1_L (R, A, DHR)
  CALL DERIV_C1_L (S, B, DHS)

! Form tensor products for R direction
  DH (1, 1) = DHR (1) * HS (1)
  DH (1, 2) = DHR (2) * HS (1)
  DH (1, 3) = DHR (1) * HS (2)
  DH (1, 4) = DHR (2) * HS (2)
  DH (1, 5) = DHR (3) * HS (1)
  DH (1, 6) = DHR (4) * HS (1)
  DH (1, 7) = DHR (3) * HS (2)
  DH (1, 8) = DHR (4) * HS (2)
  DH (1, 9) = DHR (3) * HS (3)
  DH (1, 10) = DHR (4) * HS (3)
  DH (1, 11) = DHR (3) * HS (4)
  DH (1, 12) = DHR (4) * HS (4)
  DH (1, 13) = DHR (1) * HS (3)
  DH (1, 14) = DHR (2) * HS (3)
  DH (1, 15) = DHR (1) * HS (4)
  DH (1, 16) = DHR (2) * HS (4)

! Form tensor products for S direction
  DH (2, 1) = HR (1) * DHS (1)
  DH (2, 2) = HR (2) * DHS (1)
  DH (2, 3) = HR (1) * DHS (2)
  DH (2, 4) = HR (2) * DHS (2)
  DH (2, 5) = HR (3) * DHS (1)
  DH (2, 6) = HR (4) * DHS (1)
  DH (2, 7) = HR (3) * DHS (2)
  DH (2, 8) = HR (4) * DHS (2)
  DH (2, 9) = HR (3) * DHS (3)
  DH (2, 10) = HR (4) * DHS (3)
  DH (2, 11) = HR (3) * DHS (4)
  DH (2, 12) = HR (4) * DHS (4)
  DH (2, 13) = HR (1) * DHS (3)
  DH (2, 14) = HR (2) * DHS (3)
  DH (2, 15) = HR (1) * DHS (4)
  DH (2, 16) = HR (2) * DHS (4)
END SUBROUTINE DERIV_16_R
                         
SUBROUTINE DERIV_17_Q (R, S, DH)
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-*
! LOCAL DERIVATIVES OF SERENDIPITY QUAD WITH 17 NODES
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-*
! SEE SHAPE_17_Q FOR NODE LOCATIONS
Use Precision_Module
 IMPLICIT NONE
 REAL(DP), PARAMETER   :: PT667  = 0.66666666666667d0, &
                          PT1D12 = 0.083333333333333d0, HALF = 0.5d0
 REAL(DP), INTENT(IN)  :: R, S
 REAL(DP), INTENT(OUT) :: DH (2, 17)
 REAL(DP)              :: RR, SS, RRR, SSS, RP, RM, S_P, SM, R1
  RP = 1.d0 + R     ; RM = 1.d0 - R
  S_P = 1.d0 + S     ; SM = 1.d0 - S
  RR = R * R      ; SS = S * S
  RRR = R * R * R ; SSS = S * S * S
  R1 = R - 1.
  DH (1, 1) = PT1D12 * SM * (16.d0 * RRR - 12.d0 * RR &
            - 6.d0 * R * S - 8.d0 * R + 4.d0 * SSS - S + 4.d0)
  DH (1, 2) = PT1D12 * SM * (16.d0 * RRR + 12.d0 * RR &
            - 6.d0 * R * S - 8.d0 * R - 4.d0 * SSS + S - 4.d0)
  DH (1, 3) = PT1D12 * S_P * (16.d0 * RRR + 12.d0 * RR &
            + 6.d0 * R * S - 8.d0 * R + 4.d0 * SSS - S - 4.d0)
  DH (1, 4) = PT1D12 * S_P * (16.d0 * RRR - 12.d0 * RR &
            + 6.d0 * R * S - 8.d0 * R - 4.d0 * SSS + S + 4.d0)
  DH (1, 5) = -PT667 * (1.d0 - 4.d0 * R - 3.d0 * RR + 8.d0 * RRR) * SM
  DH (1, 6) = PT667 * S * SM * S_P * ( - 1.d0 + 2.d0 * S)
  DH (1, 7) = -PT667 * (-1.d0 - 4.d0 * R + 3.d0 * RR + 8.d0 * RRR) * S_P
  DH (1, 8) = - PT667 * S * SM * S_P * (1.d0 + 2.d0 * S)
  DH (1, 9) = R * SM * (8.d0 * RR + S - 4.d0)
  DH (1, 10) = HALF * SM * S_P * (2.d0 * R - 4.d0 * SS + 1.d0)
  DH (1, 11) = R * S_P * (8.d0 * RR - S - 4.d0)
  DH (1, 12) = HALF * SM * S_P * (2.d0 * R - 1.d0 + 4.d0 * SS)
  DH (1, 13) = PT667 * (1.d0 + 4.d0 * R - 3.d0 * RR - 8.d0 * RRR) * SM
  DH (1, 14) = PT667 * S * SM * S_P * (1.d0 + 2.d0 * S)
  DH (1, 15) = -PT667 * (1.d0 - 4.d0 * R - 3.d0 * RR + 8.d0 * RRR) * S_P
  DH (1, 16) = PT667 * S * SM * S_P * (1.d0 - 2.d0 * S)
  DH (1, 17) = - 2.d0 * R * SM * S_P
  DH (2, 1) = PT1D12 * RM * (16.d0 * SSS - 12.d0 * SS &
            - 6.d0 * R * S - 8.d0 * S + 4.d0 * RRR - R + 4.d0)
  DH (2, 2) = - PT1D12 * RP * ( -16.d0 * SSS + 12.d0 * SS &
            - 6.d0 * R * S + 8.d0 * S + 4.d0 * RRR - R - 4.d0)
  DH (2, 3) = PT1D12 * RP * (16.d0 * SSS + 12.d0 * SS &
            + 6.d0 * R * S - 8.d0 * S + 4.d0 * RRR - R - 4.d0)
  DH (2, 4) = PT1D12 * R1 * ( -16.d0 * SSS - 12.d0 * SS &
            + 6.d0 * R * S + 8.d0 * S + 4.d0 * RRR - R + 4.d0)
  DH (2, 5) = PT667 * R * R1 * RP * (2.d0 * R - 1.d0)
  DH (2, 6) = - PT667 * (1.d0 - 4.d0 * S - 3.d0 * SS + 8.d0 * SSS) * RP
  DH (2, 7) = - PT667 * R * R1 * RP * (1.d0 + 2.d0 * R)
  DH (2, 8) = PT667 * ( - 1.d0 - 4.d0 * S + 3.d0 * SS + 8.d0 * SSS) * R1
  DH (2, 9) = - HALF * R1 * RP * (2.d0 * S - 1.d0 + 4.d0 * RR)
  DH (2, 10) = - S * RP * ( - 8.d0 * SS + R + 4.d0)
  DH (2, 11) = HALF * R1 * RP * ( - 2.d0 * S + 4.d0 * RR - 1.d0)
  DH (2, 12) = - S * R1 * (8.d0 * SS + R - 4.d0)
  DH (2, 13) = PT667 * R * R1 * RP * (1.d0 + 2.d0 * R)
  DH (2, 14) = -PT667 * (-1.d0 - 4.d0 * S + 3.d0 * SS + 8.d0 * SSS) * RP
  DH (2, 15) = -PT667 * R * R1 * RP * (2.d0 * R - 1.d0)
  DH (2, 16) = PT667 * (1.d0 - 4.d0 * S - 3.d0 * SS + 8.d0 * SSS) * R1
  DH (2, 17) = 2.d0 * S * R1 * RP
END SUBROUTINE DERIV_17_Q
                         
SUBROUTINE DERIV_8_20_H (R, S, T, DH, ELEM_NODES)
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *
!   LOCAL DERIVATIVES OF INTERPOLATION FUNCTIONS FOR AN
! 8 TO 20 NODE HEXAHEDRON, SEE SHAPE_8_20_H FOR TOPOLOGY FIGURE
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *
Use Precision_Module
 IMPLICIT NONE
 REAL(DP), INTENT(IN)  :: R, S, T
 REAL(DP), INTENT(OUT) :: DH (3, 20)
 INTEGER  :: ELEM_NODES (20), I1 (20), I2 (20), K, K1, K2
 REAL(DP) :: FR, FS, FT, RP, RM, S_P, SM, TP, TM, &
             RZ, SZ, TZ, DH1, DH2, DH3, FP, FM
  DATA I1 / 8 * 0, 1, 2, 3, 4, 5, 6, 7, 8, 1, 2, 3, 4 /
  DATA I2 / 8 * 0, 2, 3, 4, 1, 6, 7, 8, 5, 5, 6, 7, 8 /
  DATA FP, FM / 0.5d0, -0.5d0 /

! FOR PARAMETER DEFINITIONS SEE SUBROUTINE SHAPE_8_20_H
! R,S,T = LOCAL COORDINATES OF THE POINT   -1 <= (R,S,T) <= +1
! DH    = LOCAL DERIVATIVES OF SHAPE FUNCTIONS, 0 IF ELEM_NODES(I) = 0
! ELEM_NODES = ARRAY OF ELEMENT INCIDENCES,
!   IF ELEM_NODES(I)=0 THEN LOCAL NODE I IS NOT CONSIDERED IN ANALYSIS
! I1, I2 = CORNER NODES OF TWELVE EDGES

  RP = 0.5d0 * (1.d0 + R) ; S_P = 0.5d0 * (1.d0 + S) ; TP = 0.5d0 * (1.d0 + T)
  RM = 0.5d0 * (1.d0 - R) ; SM = 0.5d0 * (1.d0 - S) ; TM = 0.5d0 * (1.d0 - T)
  RZ = 1.d0 - R * R     ; SZ = 1.d0 - S * S     ; TZ = 1.d0 - T * T
  FR = - 2.d0 * R      ; FS = - 2.d0 * S      ; FT = - 2.d0 * T

!  DERIVATIVES OF TRI-LINEAR CORNERS
  DH (1, 1) = TP * S_P * FP
  DH (2, 1) = TP * FP * RP
  DH (3, 1) = FP * S_P * RP
  DH (1, 2) = TP * S_P * FM
  DH (2, 2) = TP * FP * RM
  DH (3, 2) = FP * S_P * RM
  DH (1, 3) = TP * SM * FM
  DH (2, 3) = TP * FM * RM
  DH (3, 3) = FP * SM * RM
  DH (1, 4) = TP * SM * FP
  DH (2, 4) = TP * FM * RP
  DH (3, 4) = FP * SM * RP
  DH (1, 5) = TM * S_P * FP
  DH (2, 5) = TM * FP * RP
  DH (3, 5) = FM * S_P * RP
  DH (1, 6) = TM * S_P * FM
  DH (2, 6) = TM * FP * RM
  DH (3, 6) = FM * S_P * RM
  DH (1, 7) = TM * SM * FM
  DH (2, 7) = TM * FM * RM
  DH (3, 7) = FM * SM * RM
  DH (1, 8) = TM * SM * FP
  DH (2, 8) = TM * FM * RP
  DH (3, 8) = FM * SM * RP

!  DERIVATIVES OF EDGE BUBBLES
  DH (1, 9) = TP * S_P * FR * 0.5d0
  DH (2, 9) = TP * FP * RZ * 0.5d0
  DH (3, 9) = FP * S_P * RZ * 0.5d0
  DH (1, 10) = TP * SZ * FM * 0.5d0
  DH (2, 10) = TP * FS * RM * 0.5d0
  DH (3, 10) = FP * SZ * RM * 0.5d0
  DH (1, 11) = TP * SM * FR * 0.5d0
  DH (2, 11) = TP * FM * RZ * 0.5d0
  DH (3, 11) = FP * SM * RZ * 0.5d0
  DH (1, 12) = TP * SZ * FP * 0.5d0
  DH (2, 12) = TP * FS * RP * 0.5d0
  DH (3, 12) = FP * SZ * RP * 0.5d0
  DH (1, 13) = TM * S_P * FR * 0.5d0
  DH (2, 13) = TM * FP * RZ * 0.5d0
  DH (3, 13) = FM * S_P * RZ * 0.5d0
  DH (1, 14) = TM * SZ * FM * 0.5d0
  DH (2, 14) = TM * FS * RM * 0.5d0
  DH (3, 14) = FM * SZ * RM * 0.5d0
  DH (1, 15) = TM * SM * FR * 0.5d0
  DH (2, 15) = TM * FM * RZ * 0.5d0
  DH (3, 15) = FM * SM * RZ * 0.5d0
  DH (1, 16) = TM * SZ * FP * 0.5d0
  DH (2, 16) = TM * FS * RP * 0.5d0
  DH (3, 16) = FM * SZ * RP * 0.5d0
  DH (1, 17) = TZ * S_P * FP * 0.5d0
  DH (2, 17) = TZ * FP * RP * 0.5d0
  DH (3, 17) = FT * S_P * RP * 0.5d0
  DH (1, 18) = TZ * S_P * FM * 0.5d0
  DH (2, 18) = TZ * FP * RM * 0.5d0
  DH (3, 18) = FT * S_P * RM * 0.5d0
  DH (1, 19) = TZ * SM * FM * 0.5d0
  DH (2, 19) = TZ * FM * RM * 0.5d0
  DH (3, 19) = FT * SM * RM * 0.5d0
  DH (1, 20) = TZ * SM * FP * 0.5d0
  DH (2, 20) = TZ * FM * RP * 0.5d0
  DH (3, 20) = FT * SM * RP * 0.5d0

!  LOOP OVER TWELVE EDGES
  DO K = 9, 20
    IF (ELEM_NODES (K)  == 0) THEN
!      ZERO EDGE BUBBLE DERIVATIVES
      DH (1, K) = 0.d0
      DH (2, K) = 0.d0
      DH (3, K) = 0.d0
    ELSE
!      ENRICH DERIVATIVES AT TWO ENDS OF EDGE
      DH1 = DH (1, K)
      DH2 = DH (2, K)
      DH3 = DH (3, K)
      K1 = I1 (K)
      K2 = I2 (K)
      DH (1, K1) = DH (1, K1) - DH1
      DH (2, K1) = DH (2, K1) - DH2
      DH (3, K1) = DH (3, K1) - DH3
      DH (1, K2) = DH (1, K2) - DH1
      DH (2, K2) = DH (2, K2) - DH2
      DH (3, K2) = DH (3, K2) - DH3
      DH (1, K) = DH1 + DH1
      DH (2, K) = DH2 + DH2
      DH (3, K) = DH3 + DH3
    END IF
  END DO
END SUBROUTINE DERIV_8_20_H
                           
SUBROUTINE DERIV_20_H (R, S, T, DH)
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *
!   LOCAL DERIVATIVES OF INTERPOLATION FUNCTIONS FOR AN
!  20 NODE HEXAHEDRON, SEE SHAPE_20_H FOR TOPOLOGY FIGURE
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *
Use Precision_Module
 IMPLICIT NONE
 REAL(DP), INTENT(IN)  :: R, S, T
 REAL(DP), INTENT(OUT) :: DH (3, 20)
 INTEGER  :: I1 (20), I2 (20), K, K1, K2
 REAL(DP) :: FR, FS, FT, RP, RM, S_P, SM, TP, TM, &
             RZ, SZ, TZ, DH1, DH2, DH3, FP, FM
  DATA I1 / 8 * 0, 1, 2, 3, 4, 5, 6, 7, 8, 1, 2, 3, 4 /
  DATA I2 / 8 * 0, 2, 3, 4, 1, 6, 7, 8, 5, 5, 6, 7, 8 /
  DATA FP, FM / 0.5d0, -0.5d0 /

! FOR PARAMETER DEFINITIONS SEE SUBROUTINE SHAPE_20_H
! R,S,T = LOCAL COORDINATES OF THE POINT   -1 <= (R,S,T) <= +1
! DH    = LOCAL DERIVATIVES OF SHAPE FUNCTIONS
! I1, I2 = CORNER NODES OF TWELVE EDGES

  RP = 0.5d0 * (1.d0 + R) ; S_P = 0.5d0 * (1.d0 + S) ; TP = 0.5d0 * (1.d0 + T)
  RM = 0.5d0 * (1.d0 - R) ; SM = 0.5d0 * (1.d0 - S) ; TM = 0.5d0 * (1.d0 - T)
  RZ = 1.d0 - R * R     ; SZ = 1.d0 - S * S     ; TZ = 1.d0 - T * T
  FR = - 2.d0 * R      ; FS = - 2.d0 * S      ; FT = - 2.d0 * T

!  DERIVATIVES OF TRI-LINEAR CORNERS
  DH (1, 1) = TP * S_P * FP
  DH (2, 1) = TP * FP * RP
  DH (3, 1) = FP * S_P * RP
  DH (1, 2) = TP * S_P * FM
  DH (2, 2) = TP * FP * RM
  DH (3, 2) = FP * S_P * RM
  DH (1, 3) = TP * SM * FM
  DH (2, 3) = TP * FM * RM
  DH (3, 3) = FP * SM * RM
  DH (1, 4) = TP * SM * FP
  DH (2, 4) = TP * FM * RP
  DH (3, 4) = FP * SM * RP
  DH (1, 5) = TM * S_P * FP
  DH (2, 5) = TM * FP * RP
  DH (3, 5) = FM * S_P * RP
  DH (1, 6) = TM * S_P * FM
  DH (2, 6) = TM * FP * RM
  DH (3, 6) = FM * S_P * RM
  DH (1, 7) = TM * SM * FM
  DH (2, 7) = TM * FM * RM
  DH (3, 7) = FM * SM * RM
  DH (1, 8) = TM * SM * FP
  DH (2, 8) = TM * FM * RP
  DH (3, 8) = FM * SM * RP

!  DERIVATIVES OF EDGE BUBBLES
  DH (1, 9) = TP * S_P * FR * 0.5d0
  DH (2, 9) = TP * FP * RZ * 0.5d0
  DH (3, 9) = FP * S_P * RZ * 0.5d0
  DH (1, 10) = TP * SZ * FM * 0.5d0
  DH (2, 10) = TP * FS * RM * 0.5d0
  DH (3, 10) = FP * SZ * RM * 0.5d0
  DH (1, 11) = TP * SM * FR * 0.5d0
  DH (2, 11) = TP * FM * RZ * 0.5d0
  DH (3, 11) = FP * SM * RZ * 0.5d0
  DH (1, 12) = TP * SZ * FP * 0.5d0
  DH (2, 12) = TP * FS * RP * 0.5d0
  DH (3, 12) = FP * SZ * RP * 0.5d0
  DH (1, 13) = TM * S_P * FR * 0.5d0
  DH (2, 13) = TM * FP * RZ * 0.5d0
  DH (3, 13) = FM * S_P * RZ * 0.5d0
  DH (1, 14) = TM * SZ * FM * 0.5d0
  DH (2, 14) = TM * FS * RM * 0.5d0
  DH (3, 14) = FM * SZ * RM * 0.5d0
  DH (1, 15) = TM * SM * FR * 0.5d0
  DH (2, 15) = TM * FM * RZ * 0.5d0
  DH (3, 15) = FM * SM * RZ * 0.5d0
  DH (1, 16) = TM * SZ * FP * 0.5d0
  DH (2, 16) = TM * FS * RP * 0.5d0
  DH (3, 16) = FM * SZ * RP * 0.5d0
  DH (1, 17) = TZ * S_P * FP * 0.5d0
  DH (2, 17) = TZ * FP * RP * 0.5d0
  DH (3, 17) = FT * S_P * RP * 0.5d0
  DH (1, 18) = TZ * S_P * FM * 0.5d0
  DH (2, 18) = TZ * FP * RM * 0.5d0
  DH (3, 18) = FT * S_P * RM * 0.5d0
  DH (1, 19) = TZ * SM * FM * 0.5d0
  DH (2, 19) = TZ * FM * RM * 0.5d0
  DH (3, 19) = FT * SM * RM * 0.5d0
  DH (1, 20) = TZ * SM * FP * 0.5d0
  DH (2, 20) = TZ * FM * RP * 0.5d0
  DH (3, 20) = FT * SM * RP * 0.5d0

!  LOOP OVER TWELVE EDGES
  DO K = 9, 20
!    ENRICH DERIVATIVES AT TWO ENDS OF EDGE
    DH1 = DH (1, K)
    DH2 = DH (2, K)
    DH3 = DH (3, K)
    K1 = I1 (K)
    K2 = I2 (K)
    DH (1, K1) = DH (1, K1) - DH1
    DH (2, K1) = DH (2, K1) - DH2
    DH (3, K1) = DH (3, K1) - DH3
    DH (1, K2) = DH (1, K2) - DH1
    DH (2, K2) = DH (2, K2) - DH2
    DH (3, K2) = DH (3, K2) - DH3
    DH (1, K) = DH1 + DH1
    DH (2, K) = DH2 + DH2
    DH (3, K) = DH3 + DH3
  END DO
END SUBROUTINE DERIV_20_H
                           
SUBROUTINE DERIV_2_L (R, DH)
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *
! DERIVATIVES OF A 2 NODE LINE ELEMENT
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *
Use Precision_Module
 IMPLICIT NONE
 REAL(DP), INTENT(IN)  :: R
 REAL(DP), INTENT(OUT) :: DH (2)

! R IS UNIT COORD.   R=-1  1------------2  R=1 (unused)

  DH (1) = - 0.5d0
  DH (2) = 0.5d0
END SUBROUTINE DERIV_2_L
                        
SUBROUTINE DERIV_3_T (S, T, DH)
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-*
! LOCAL DERIVATIVES OF A THREE NODE UNIT TRIANGLE
!            SEE SUBROUTINE SHAPE_3_T
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-*
Use Precision_Module
 IMPLICIT NONE
 REAL(DP), INTENT(IN)  :: S, T
 REAL(DP), INTENT(OUT) :: DH (2, 3)

! S,T     = LOCAL COORDINATES OF THE POINT
! DH(1,K) = DH(K)/DS
! DH(2,K) = DH(K)/DT
! NODAL COORDS ARE :  1-(0,0) 2-(1,0) 3-(0,1)

  DH (1, 1) = - 1.d0
  DH (1, 2) = 1.d0
  DH (1, 3) = 0.d0
  DH (2, 1) = - 1.d0
  DH (2, 2) = 0.d0
  DH (2, 3) = 1.d0
END SUBROUTINE DERIV_3_T
                        
!SUBROUTINE DERIV_4_12_Q (R, S, DH, ELEM_NODES)
!! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-*
!!    DERIVATIVES OF 4 TO 12 NODE QUADRILATERAL
!! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-*
!Use Precision_Module
! IMPLICIT NONE
! REAL(DP), INTENT(IN)  :: R, S
! INTEGER,  INTENT(IN)  :: ELEM_NODES (12)
! REAL(DP), INTENT(OUT) :: DH (2, 12)
!  DH = 0.d0
!  STOP 'CONSTRUCT DERIV_4_12_Q FROM SHAPE_4_12_Q'
!END SUBROUTINE DERIV_4_12_Q
                           
SUBROUTINE DERIV_4_L (X, DH)
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-*
! CALCULATE DERIVATIVE OF A 4 NODE LINE ELEMENT
!         IN NATURAL COORDINATES
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-*
Use Precision_Module
 IMPLICIT NONE
 REAL(DP), INTENT(IN)  :: X
 REAL(DP), INTENT(OUT) :: DH (4)

! DH = LOCAL DERIVATIVE OF SHAPE_4_L
! X = LOCAL COORDIATE OF POINT, -1 <= X <= +1
! LOCAL NODE COORD. ARE    1----2---3----4

  DH (1) = (1.d0 + 18.d0 * X - 27.d0 * X * X) / 16.d0
  DH (2) = 9.d0 / 16.d0 * ( - 3.d0 - 2.d0 * X + 9.d0 * X * X)
  DH (3) = 9.d0 / 16.d0 * (3.d0 - 2.d0 * X - 9.d0 * X * X)
  DH (4) = ( - 1.d0 + 18.d0 * X + 27.d0 * X * X) / 16.d0
END SUBROUTINE DERIV_4_L
                        
SUBROUTINE DERIV2_4_L (X, D2H)
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-*
! CALCULATE SECOND DERIVATIVE OF A 4 NODE LINE ELEMENT
!         IN NATURAL COORDINATES
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-*
Use Precision_Module
 IMPLICIT NONE
 REAL(DP), INTENT(IN)  :: X
 REAL(DP), INTENT(OUT) :: D2H (4)

! D2H = SECOND LOCAL DERIVATIVE OF SHAPE_4_L
! X = LOCAL COORDINATE OF POINT, -1 <= X <= +1
! LOCAL NODE COORD. ARE    1----2---3----4

  D2H (1) = (18.d0 - 54.d0 * X) / 16.d0
  D2H (2) = 9.d0 / 16.d0 * ( - 2.d0 + 18.d0 * X)
  D2H (3) = 9.d0 / 16.d0 * ( - 2.d0 - 18.d0 * X)
  D2H (4) = (18.d0 + 54.d0 * X) / 16.d0
END SUBROUTINE DERIV2_4_L
                         
SUBROUTINE DERIV_4_P (R, S, T, DH)
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-*
! LOCAL DERIVATIVE OF A 4 NODE TETRAHEDRON (PYRAMID)
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-*
Use Precision_Module
 IMPLICIT NONE
 REAL(DP), INTENT(IN)  :: R, S, T
 REAL(DP), INTENT(OUT) :: DH (3, 4)

  DH (1, 1) = - 1.d0
  DH (1, 2) = 1.d0
  DH (1, 3) = 0.d0
  DH (1, 4) = 0.d0
  DH (2, 1) = - 1.d0
  DH (2, 2) = 0.d0
  DH (2, 3) = 1.d0
  DH (2, 4) = 0.d0
  DH (3, 1) = - 1.d0
  DH (3, 2) = 0.d0
  DH (3, 3) = 0.d0
  DH (3, 4) = 1.d0
END SUBROUTINE DERIV_4_P
                        
SUBROUTINE DERIV_4_Q (R, S, DELTA)
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-*
! LOCAL DERIVATIVES OF THE SHAPE FUNCTIONS FOR AN
! PARAMETRIC QUADRILATERAL WITH FOUR NODES
!            SEE SHAPE_4_Q
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-*
Use Precision_Module
 IMPLICIT NONE
 REAL(DP), INTENT(IN)  :: R, S
 REAL(DP), INTENT(OUT) :: DELTA (2, 4)
 REAL(DP)              :: RP, RM, S_P, SM

! DELTA(1,I) = DH/DR
! DELTA(2,I) = DH/DS
! H          = LOCAL INTERPOLATION FUNCTIONS
! (R,S)      = A POINT IN THE LOCAL COORDINATES
! HERE D(H(I))/DR = 0.25d0*R(I)*(1+S*S(I)), ETC.

  RP = 1.d0 + R ; RM = 1.d0 - R
  S_P = 1.d0 + S ; SM = 1.d0 - S
  DELTA (1, 1) = - 0.25d0 * SM
  DELTA (1, 2) = 0.25d0 * SM
  DELTA (1, 3) = 0.25d0 * S_P
  DELTA (1, 4) = - 0.25d0 * S_P
  DELTA (2, 1) = - 0.25d0 * RM
  DELTA (2, 2) = - 0.25d0 * RP
  DELTA (2, 3) = 0.25d0 * RP
  DELTA (2, 4) = 0.25d0 * RM
END SUBROUTINE DERIV_4_Q
                        
SUBROUTINE DERIV2_4_Q (R, S, D2L_H)
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-*
! LOCAL SECOND DERIVATIVES OF THE SHAPE FUNCTIONS FOR 
!    A PARAMETRIC QUADRILATERAL WITH FOUR NODES
!                   H,RR H,RS H,SS
!            SEE SHAPE_4_Q, DERIV_4_Q
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-*
Use Precision_Module
 IMPLICIT NONE
 REAL(DP), INTENT(IN)  :: R, S
 REAL(DP), INTENT(OUT) :: D2L_H (3, 4)

! D2L_H(1,I) = D2 H/DR DR FOR NODE I
! D2L_H(2,I) = D2 H/DR DS
! D2L_H(3,I) = D2 H/DS DS
! H          = LOCAL INTERPOLATION FUNCTIONS
! (R,S)      = A POINT IN THE LOCAL COORDINATES
! HERE D(H(I))/DR = 0.25d0*R(I)*(1+S*S(I)), ETC.

  D2L_H (1, 1:4) = 0.d0 
  D2L_H (2, 1:4) = (/ -0.25d0, 0.25d0, 0.25d0, -0.25d0 /)
  D2L_H (3, 1:4) = 0.d0 
END SUBROUTINE DERIV2_4_Q
                        
SUBROUTINE DERIV_4_R (R, S, A, B, DH)
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *
!       FIRST DERIVATIVES OF A
! C0 RECTANGULAR ELEMENT IN UNIT COORDINATES
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *
Use Precision_Module
 IMPLICIT NONE
 REAL(DP), INTENT(IN)  :: R, S, A, B
 REAL(DP), INTENT(OUT) :: DH (2, 4)

! DOF ARE W AT EACH NODE (N_G_DOF=1)
! X // R, Y // S.              S
! A = PHYSICAL LENGTH IN X     4 -------- 3
! B = PHYSICAL LENGTH IN Y     I          I
! R,S = LOCAL UNIT COORDS      I          I
! 1@(0,0), 3@(1,1)             1 -------- 2 ->R

  DH (1, 1) = (S - 1.d0) / A
  DH (1, 2) = (1.d0 - S) / A
  DH (1, 3) = S / A
  DH (1, 4) = - S / A
  DH (2, 1) = (R - 1.d0) / B
  DH (2, 2) = - R / B
  DH (2, 3) = R / B
  DH (2, 4) = (1.d0 - R) / B
END SUBROUTINE DERIV_4_R
                        
SUBROUTINE DERIV_4_T (R, S, DH)
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-*
! LOCAL DERIVATIVES OF A
! UNIT TRIANGLE WITH CENTROID (4TH) NODE
!   This is untested
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-*
Use Precision_Module
 IMPLICIT NONE
 REAL(DP), INTENT(IN)  :: R, S
 REAL(DP), INTENT(OUT) :: DH (2, 4)

! R,S = UNIT COORDINATES, SEE SHAPE_4_T
! DH = DERIVATIVES OF INTERPOLATION FUNCTIONS
! NODE 1@(0,0), 2@(1,0), 3@(0,1), 4@(1/3,1/3)

  DH (1, 4) = 27.d0 * S * (1.d0 - R - S) - 27.d0 * R * S
  DH (2, 4) = 27.d0 * R * (1.d0 - R - S) - 27.d0 * R * S
  DH (1, 1) = - 1.d0 - DH (1, 4) / 3.d0
  DH (2, 1) = - 1.d0 - DH (2, 4) / 3.d0
  DH (1, 2) = 1.d0 - DH (1, 4) / 3.d0
  DH (2, 2) = - DH (2, 4) / 3.d0
  DH (1, 3) = - DH (1, 4) / 3.d0
  DH (2, 3) = 1.d0 - DH (2, 4) / 3.d0
END SUBROUTINE DERIV_4_T
                        
SUBROUTINE DERIV_7_T (S, T, DH)
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-*
! LOCAL DERIVATIVES FOR A SIX NODE UNIT TRIANGLE
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-*
Use Precision_Module
 IMPLICIT NONE
 REAL(DP), INTENT(IN)  :: S, T
 REAL(DP), INTENT(OUT) :: DH (2, 7)

! S,T = LOCAL COORDINATES, SEE SHAPE_7_T
! DH = LOCAL DERIVATIVES OF SHAPE FUNCTIONS
! DH(1,K) = DH(K)/DS,   DH(2,K)=DH(K)/DT

  DH = 0.d0
  STOP 'CONSTRUCT DERIV_7_T FROM SHAPE_7_T'
END SUBROUTINE DERIV_7_T
                        
SUBROUTINE DERIV_8_H (R, S, T, DH)
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-*
! LOCAL DERIVATIVES FOR EIGHT NODE HEXAHEDRON
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-*
Use Precision_Module
 IMPLICIT NONE
 REAL(DP), INTENT(IN)  :: R, S, T
 REAL(DP), INTENT(OUT) :: DH (3, 8)
 REAL(DP)              :: RP, RM, S_P, SM, TP, TM

! R,S,T = LOCAL COORDINATES OF THE POINT
! DH(1,K)=DH/DR, DH(2,K)=DH/DS, DH(3,K)=DH/DT
! H = ELEMENT SHAPE FUNCTIONS, SEE SHAPE_8_H

  RP = 1.d0 + R ; RM = 1.d0 - R
  S_P = 1.d0 + S ; SM = 1.d0 - S
  TP = 1.d0 + T ; TM = 1.d0 - T
  DH (1, 1) = 0.125d0 * S_P * TP
  DH (2, 1) = 0.125d0 * RP * TP
  DH (3, 1) = 0.125d0 * RP * S_P
  DH (1, 2) = 0.125d0 * SM * TP
  DH (2, 2) = - 0.125d0 * RP * TP
  DH (3, 2) = 0.125d0 * RP * SM
  DH (1, 3) = 0.125d0 * SM * TM
  DH (2, 3) = - 0.125d0 * RP * TM
  DH (3, 3) = - 0.125d0 * RP * SM
  DH (1, 4) = 0.125d0 * S_P * TM
  DH (2, 4) = 0.125d0 * RP * TM
  DH (3, 4) = - 0.125d0 * RP * S_P
  DH (1, 5) = - 0.125d0 * S_P * TP
  DH (2, 5) = 0.125d0 * RM * TP
  DH (3, 5) = 0.125d0 * RM * S_P
  DH (1, 6) = - 0.125d0 * SM * TP
  DH (2, 6) = - 0.125d0 * RM * TP
  DH (3, 6) = 0.125d0 * RM * SM
  DH (1, 7) = - 0.125d0 * SM * TM
  DH (2, 7) = - 0.125d0 * RM * TM
  DH (3, 7) = - 0.125d0 * RM * SM
  DH (1, 8) = - 0.125d0 * S_P * TM
  DH (2, 8) = 0.125d0 * RM * TM
  DH (3, 8) = - 0.125d0 * RM * S_P
END SUBROUTINE DERIV_8_H
                        
SUBROUTINE DERIV_9_Q (R, S, DH)
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *
! LOCAL DERIVATIVES FOR 9-NODED QUAD
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *
! SEE SHAPE_9_Q FOR TOPOLOGY
Use Precision_Module
 IMPLICIT NONE
 REAL(DP), INTENT(IN)  :: R, S
 REAL(DP), INTENT(OUT) :: DH (2, 9)
 REAL(DP)              :: RP, RM, S_P, SM, R2P1, R2M1, S2P1, S2M1

  RM = R - 1.D0        ; SM = S - 1.D0
  RP = R + 1.D0        ; S_P = S + 1.D0
  S2P1 = S + S + 1.D0  ; S2M1 = S + S - 1.D0
  R2P1 = R + R + 1.D0  ; R2M1 = R + R - 1.D0
  DH (1, 1) = 0.25D0 * S * SM * R2M1
  DH (1, 2) = 0.25D0 * S * SM * R2P1
  DH (1, 3) = 0.25D0 * S * S_P * R2P1
  DH (1, 4) = 0.25D0 * S * S_P * R2M1
  DH (1, 5) = - S * SM * R
  DH (1, 6) = - 0.5D0 * S_P * SM * R2P1
  DH (1, 7) = - S * S_P * R
  DH (1, 8) = - 0.5D0 * S_P * SM * R2M1
  DH (1, 9) = 2.D0 * S_P * SM * R
  DH (2, 1) = 0.25D0 * S2M1 * R * RM
  DH (2, 2) = 0.25D0 * S2M1 * R * RP
  DH (2, 3) = 0.25D0 * S2P1 * R * RP
  DH (2, 4) = 0.25D0 * S2P1 * R * RM
  DH (2, 5) = - 0.5D0 * S2M1 * RP * RM
  DH (2, 6) = - S * R * RP
  DH (2, 7) = - 0.5D0 * S2P1 * RP * RM
  DH (2, 8) = - S * R * RM
  DH (2, 9) = 2.D0 * S * RP * RM
END SUBROUTINE DERIV_9_Q

SUBROUTINE DERIV2_9_Q (R, S, D2L_H)
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *
! LOCAL SECOND DERIVATIVES FOR 9-NODED QUAD
!                  H,RR H,RS H,SS
!             SEE SHAPE_9_Q FOR TOPOLOGY
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *
Use Precision_Module
 IMPLICIT NONE
 REAL(DP), INTENT(IN)  :: R, S
 REAL(DP), INTENT(OUT) :: D2L_H (3, 9)
 REAL(DP)              :: RP, RM, S_P, SM, R2P1, R2M1, S2P1, S2M1

! D2L_H(1,I) = D2 H/DR DR FOR NODE I
! D2L_H(2,I) = D2 H/DR DS
! D2L_H(3,I) = D2 H/DS DS
! H          = LOCAL INTERPOLATION FUNCTIONS
! (R,S)      = A POINT IN THE LOCAL COORDINATES

  RM = R - 1.D0        ; SM = S - 1.D0
  RP = R + 1.D0        ; S_P = S + 1.D0
  S2P1 = S + S + 1.D0  ; S2M1 = S + S - 1.D0
  R2P1 = R + R + 1.D0  ; R2M1 = R + R - 1.D0

  D2L_H (1, 1) = 0.5D0 * S * SM
  D2L_H (1, 2) = 0.5D0 * S * SM 
  D2L_H (1, 3) = 0.5D0 * S * S_P
  D2L_H (1, 4) = 0.5D0 * S * S_P 
  D2L_H (1, 5) = - S * SM 
  D2L_H (1, 6) = - S_P * SM 
  D2L_H (1, 7) = - S * S_P 
  D2L_H (1, 8) = - S_P * SM 
  D2L_H (1, 9) = 2.D0 * S_P * SM 

  D2L_H (2, 1) = 0.25D0 * S2M1 * R2M1   
  D2L_H (2, 2) = 0.25D0 * S2M1 * R2P1   
  D2L_H (2, 3) = 0.25D0 * S2P1 * R2P1   
  D2L_H (2, 4) = 0.25D0 * S2P1 * R2M1   
  D2L_H (2, 5) = - S2M1 * R
  D2L_H (2, 6) = - S * R2P1    
  D2L_H (2, 7) = - S2P1 * R
  D2L_H (2, 8) = - S * R2M1   
  D2L_H (2, 9) = 4.D0 * S * R

  D2L_H (3, 1) = 0.5D0 * R * RM
  D2L_H (3, 2) = 0.5D0 * R * RP
  D2L_H (3, 3) = 0.5D0 * R * RP
  D2L_H (3, 4) = 0.5D0 * R * RM
  D2L_H (3, 5) = - RP * RM
  D2L_H (3, 6) = - R * RP
  D2L_H (3, 7) = - RP * RM
  D2L_H (3, 8) = - R * RM
  D2L_H (3, 9) = 2.D0 * RP * RM
END SUBROUTINE DERIV2_9_Q
                        
SUBROUTINE DERIV_C0_L (R, A, DH)
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-*
! FIRST DERIVATIVE OF C0 LINE ELEMENT IN UNIT COORDINATES
!                 (SEE SHAPE_C0_L)
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-*
Use Precision_Module
 IMPLICIT NONE
 REAL(DP), INTENT(IN)  :: R, A
 REAL(DP), INTENT(OUT) :: DH (2)

! A = PHYSICAL LENGTH OF ELEMENT       1 -------- 2 --> R
! R = LOCAL COORDINATE OF POINT        R=0        R=1
! DH = FIRST PHYSICAL DERIVATIVES OF H
! D()/DX = D()/DR DR/DX = 1/A * D()/DR
! DOF ARE NODAL VALUES ONLY

  DH (1) = - 1.d0 / A
  DH (2) = 1.d0 / A
END SUBROUTINE DERIV_C0_L
                         
SUBROUTINE DERIV_C1_L (R, A, DH)
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-*
! FIRST DERIVATIVES OF CUBIC HERMITE IN UNIT COORDINATES
!                 (A C1 ELEMENT, SEE SHAPE_C1_L)
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-*
Use Precision_Module
 IMPLICIT NONE
 REAL(DP), INTENT(IN)  :: R, A
 REAL(DP), INTENT(OUT) :: DH (4)

! A = PHYSICAL LENGTH OF ELEMENT       1 -------- 2 --> R
! R = LOCAL COORDINATE OF POINT        R=0        R=1
! DH = FIRST PHYSICAL DERIVATIVES OF H
! D()/DX = D()/DR DR/DX = 1/A * D()/DR
! DOF ARE FUNCTION AND SLOPE, D/DX, AT EACH NODE

  DH (1) = 6.d0 * (R * R - R) / A
  DH (2) = 1.d0 - 4.d0 * R + 3.d0 * R * R
  DH (3) = 6.d0 * (R - R * R) / A
  DH (4) = 3.d0 * R * R - 2.d0 * R
END SUBROUTINE DERIV_C1_L
                         
SUBROUTINE DERIV_C1_L2_NC (R, A, DH)
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-*
! FIRST DERIVATIVES OF CUBIC HERMITE IN NATURAL COORDINATES
!                 (A C1 ELEMENT, SEE SHAPE_C1_L2_NC)
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-*
Use Precision_Module
 IMPLICIT NONE
 REAL(DP), INTENT(IN)  :: R, A
 REAL(DP), INTENT(OUT) :: DH (4)

! A = PHYSICAL LENGTH OF ELEMENT       1 ----+---- 2 --> R
! R = LOCAL COORDINATE OF POINT        R=-1   0   R=1
! DH = FIRST PHYSICAL DERIVATIVES OF H
! D()/DX = D()/DR DR/DX = 2/A * D()/DR
! DOF ARE FUNCTION AND SLOPE, D/DX, AT EACH NODE

  DH(1) = (-3 + 3*R*R) * 0.5d0 / A
  DH(2) = (-1 - 2*R + 3*R*R) * 0.25d0
  DH(3) = ( 3 - 3*R*R) * 0.5d0 / A
  DH(4) = (-1 + 2*R + 3*R*R) * 0.25d0
END SUBROUTINE DERIV_C1_L2_NC
                         
SUBROUTINE DERIV2_C1_L (R, A, D2H)
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-*
! 2ND DERIVATIVES OF CUBIC HERMITE IN UNIT COORDINATES
!          (A C1 ELEMENT, SEE SHAPE_C1L & DERIV_C1_L)
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-*
Use Precision_Module
 IMPLICIT NONE
 REAL(DP), INTENT(IN)  :: R, A
 REAL(DP), INTENT(OUT) :: D2H (4)

! A = PHYSICAL LENGTH OF ELEMENT       1 -------- 2 --> R
! R = LOCAL COORDINATE OF POINT        R=0        R=1
! D2H = SECOND DERIVATIVES OF H
! D()/DX = D()/DR DR/DX = 1/A * D()/DR
! DOF ARE VALUE, SLOPE AT EACH END (WRT X)

  D2H (1) = 6.d0 * (R + R - 1.d0) / A**2
  D2H (2) = ( - 4.d0 + 6.d0 * R) / A
  D2H (3) = 6.d0 * (1.d0 - R - R) / A**2
  D2H (4) = (6.d0 * R - 2.d0) / A
END SUBROUTINE DERIV2_C1_L
                          
SUBROUTINE DERIV2_C1_L2_NC (R, A, D2H)
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-*
! 2ND DERIVATIVES OF CUBIC HERMITE IN UNIT COORDINATES
!   (A C1 ELEMENT, SEE SHAPE_C1_L2_NC & DERIV_C1_L2_NC)
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-*
Use Precision_Module
 IMPLICIT NONE
 REAL(DP), INTENT(IN)  :: R, A
 REAL(DP), INTENT(OUT) :: D2H (4)

! A = PHYSICAL LENGTH OF ELEMENT       1 ----+---- 2 --> R
! R = LOCAL COORDINATE OF POINT        R=-1   0   R=1
! D2H = SECOND DERIVATIVES OF H
! D()/DX = D()/DR DR/DX = 1/A * D()/DR
! DOF ARE VALUE, SLOPE AT EACH END (WRT X)

  D2H(1) = ( 6*R) / (A*A)
  D2H(2) = (-2 + 6*R) * 0.5d0 / A
  D2H(3) = (-6*R) / (A*A)
  D2H(4) = ( 2 + 6*R) * 0.5d0 / A
END SUBROUTINE DERIV2_C1_L2_NC
                          
SUBROUTINE DERIV_C2_L (R, A, DH)
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *
! FIRST DERIVATIVES OF FIFTH ORDER HERMITE ELEMENT IN UNIT
!    COORDINATES ( A C2 ELEMENT, SEE SURR SHAPE_C2_L )
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *
Use Precision_Module
 IMPLICIT NONE
 REAL(DP), INTENT(IN)  :: R, A
 REAL(DP), INTENT(OUT) :: DH (6)
 REAL(DP)              :: R3

! A = PHYSICAL LENGTH OF ELEMENT       1 -------- 2 --> R
! R = LOCAL COORDINATE OF POINT        R=0        R=1
! DH = FIRST PHYSICAL DERIVATIVES OF H
! D()/DX = D()/DR DR/DX = 1/A * D()/DR
! DOF ARE VALUE, SLOPE, CURVATURE AT EACH NODE (WRT X)

  R3 = R * R * R
  DH (1) = 30.d0 * (2.d0 * R3 - R * R - R3 * R) / A
  DH (2) = 1.d0 - 18.d0 * R * R + 32.d0 * R3 - 15.d0 * R3 * R
  DH (3) = 0.5d0*(2.d0 * R - 9.d0 * R * R + 12.d0 * R3 - 5.d0*R3*R)*A
  DH (4) = 30.d0 * R * R * (1.d0 - 2.d0 * R + R * R) / A
  DH (5) = 28.d0 * R3 - 15.d0 * R3 * R - 12.d0 * R * R
  DH (6) = 0.5d0 * R * R * (3.d0 - 8.d0 * R + 5.d0 * R * R) * A
END SUBROUTINE DERIV_C2_L
                         
SUBROUTINE DERIV2_C2_L (R, A, D2H)
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *
! 2ND DERIVATIVES OF FIFTH ORDER HERMITE ELEMENT IN UNIT
! COORDINATES ( A C2 ELEMENT, SEE SUBR SHAPE_C2_L & DERIV_C2_L )
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *
Use Precision_Module
 IMPLICIT NONE
 REAL(DP), INTENT(IN)  :: R, A
 REAL(DP), INTENT(OUT) :: D2H (6)
 REAL(DP)              :: R3

! A = PHYSICAL LENGTH OF ELEMENT       1 -------- 2 --> R
! R = LOCAL COORDINATE OF POINT        R=0        R=1
! D2H = SECOND DERIVATIVES OF H
! D()/DX = D()/DR DR/DX = 1/A * D()/DR
! DOF ARE VALUE, SLOPE, 2ND DERIV AT EACH END (WRT X)

  R3 = R * R * R
  D2H (1) = 30.d0 * (6.d0 * R**2 - 2.d0 * R - 4.d0 * R3) / A**2
  D2H (2) = ( - 36.d0 * R + 96.d0 * R * R - 60.d0 * R3) / A
  D2H (3) = 0.5d0 * (2.d0 - 18.d0 * R + 36.d0 * R * R - 20.d0 * R3)
  D2H (4) = 30.d0 * (2.d0 * R - 6.d0 * R * R + 4.d0 * R3) / A**2
  D2H (5) = (84.d0 * R * R - 60.d0 * R3 - 24.d0 * R) / A
  D2H (6) = 0.5d0 * (6.d0 * R - 24.d0 * R * R + 20.d0 * R3)
END SUBROUTINE DERIV2_C2_L
                          
SUBROUTINE DERIV3_C2_L (R, A, D3H)
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *
! 3RD DERIVATIVES OF FIFTH ORDER HERMITE ELEMENT IN UNIT
! COORDINATES (A C2 ELEMENT, SEE SUBR SHAPE_C2_L & DERIV_C2_L)
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *
Use Precision_Module
 IMPLICIT NONE
 REAL(DP), INTENT(IN)  :: R, A
 REAL(DP), INTENT(OUT) :: D3H (6)

! A = PHYSICAL LENGTH OF ELEMENT       1 -------- 2 --> R
! R = LOCAL COORDINATE OF POINT        R=0        R=1
! D3H = THIRD DERIVATIVES OF H
! D()/DX = D()/DR DR/DX = 1/A * D()/DR
! DOF ARE VALUE, SLOPE, 2ND DERIV AT EACH END (WRT X)

 !D3H (1) = 30.d0 * (6.d0 * R**2 - 2.d0 * R - 4.d0 * R3) / A**2
  D3H (1) = 30.d0 * (12.d0 * R - 2.d0 - 12.d0 * R * R) / A**3  !b
 !D3H (2) = ( - 36.d0 * R + 96.d0 * R * R - 60.d0 * R3) / A
  D3H (2) = ( - 36.d0 + 192.d0 * R - 180.d0 * R * R) / A**2  !b
 !D3H (3) = 0.5d0 * (2.d0 - 18.d0 * R + 36.d0 * R * R - 20.d0 * R3)
  D3H (3) = 0.5d0 * (-18.d0 + 72.d0 * R - 60.d0 * R * R) / A  !b
 !D3H (4) = 30.d0 * (2.d0 * R - 6.d0 * R * R + 4.d0 * R3) / A**2
  D3H (4) = 30.d0 * (2.d0 - 12.d0 * R + 12.d0 * R * R) / A**3  !b
 !D3H (5) = (84.d0 * R * R - 60.d0 * R3 - 24.d0 * R) / A
  D3H (5) = (168.d0 * R - 180.d0 * R * R - 24.d0) / A**2  !b
 !D3H (6) = 0.5d0 * (6.d0 * R - 24.d0 * R * R + 20.d0 * R3)
  D3H (6) = 0.5d0 * (6.d0 - 48.d0 * R + 60.d0 * R * R) / A  !b
END SUBROUTINE DERIV3_C2_L
                          
SUBROUTINE DERIV_C3_L (R, A, DH)
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-*
! FIRST DERIVATIVES FOR 7TH ORDER HERMITE LINE ELEMENT
!            ( A C3 ELEMENT, SEE SHAPE_C3_L )
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-*
Use Precision_Module
 IMPLICIT NONE
 REAL(DP), INTENT(IN)  :: R, A
 REAL(DP), INTENT(OUT) :: DH (8)
 REAL(DP)              :: R2, R3, R4

! A = PHYSICAL LENGTH OF ELEMENT       1 -------- 2 --> R
! R = LOCAL COORDINATE OF POINT        R=0        R=1
! DH = FIRST PHYSICAL DERIVATIVES OF H
! D()/DX = D()/DR DR/DX = 1/A * D()/DR
! DOF ARE VALUE, SLOPE, 2ND, 3RD DERIV AT EACH END (WRT X)

  R2 = R * R
  R3 = R2 * R
  R4 = R2 * R2
  DH (1) = 140.d0 * R3 * (R3 - 3.d0 * R2 + 3.d0 * R - 1.d0) / A
  DH (2) = 1.d0 - R3 * (80.d0 - 225.d0 * R + 216.d0 * R2 - 70.d0 * R3)
  DH (3) = R * (1.d0 - 20.d0*R2 + 50.d0*R3 - 45.d0*R4 + 14.d0*R2*R3) * A
  DH (4) = R2 * (3.d0 - 16.d0 * R + 30.d0 * R2 - 24.d0 * R3 &
         + 7.d0 * R4) * A * A / 6.
  DH (5) = 140.d0 * R3 * (1.d0 - 3.d0 * R + 3.d0 * R2 - R3) / A
  DH (6) = R3 * (70.d0 * R3 - 204.d0 * R2 + 195.d0 * R - 60.)
  DH (7) = R3 * (10.d0 - 35.d0 * R + 39.d0 * R2 - 14.d0 * R3) * A
  DH (8) = R3 * (7.d0 * R3 - 18.d0 * R2 + 15.d0 * R - 4.d0) * A * A / 6.
END SUBROUTINE DERIV_C3_L
                         
SUBROUTINE DERIV2_C3_L (R, A, D2H)
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-*
! 2ND DERIV FOR SEVENTH ORDER HERMITE LINE ELEMENT
!      ( A C3 ELEMENT, SEE SHAPE_C3_L & DERIV_C3_L )
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-*
Use Precision_Module
 IMPLICIT NONE
 REAL(DP), INTENT(IN)  :: R, A
 REAL(DP), INTENT(OUT) :: D2H (8)
 REAL(DP)              :: R2, R3, R4

! A = PHYSICAL LENGTH OF ELEMENT       1 -------- 2 --> R
! R = LOCAL COORDINATE OF POINT        R=0        R=1
! D()/DX = D()/DR DR/DX = 1/A * D()/DR
! D2H = SECOND DERIVATIVES OF SHAPE FUNCTIONS ARRAY
! DOF ARE VALUE, SLOPE, 2ND, 3RD DERIV AT EACH END (WRT X)

  R2 = R*R
  R3 = R2*R
  R4 = R2*R2
  D2H (1) = 140.d0*(6.d0*R2*R3 - 15.d0*R4 + 12.d0*R3 - 3.d0*R2) / A**2
  D2H (2) = - (240.d0*R2 - 900.d0*R3 + 1080.d0*R4 - 420.d0*R3*R2) / A
  D2H (3) = (1.d0 - 60.d0*R2 + 200.d0*R3 - 225.d0*R4 + 84.d0*R2*R3)
  D2H (4) = (R - 8.d0*R2 + 20.d0*R3 - 20.d0*R4 + 7.d0*R2*R3)*A
  D2H (5) = 140.d0*(3.d0*R2 - 12.d0*R3 + 15.d0*R4 - 6.d0*R2*R3) / A**2
  D2H (6) = (420.d0*R3*R2 - 1020.d0*R4 + 780.d0*R3 - 180.d0*R2) / A
  D2H (7) = (30.d0*R2 - 140.d0*R3 + 195.d0*R4 - 84.d0*R3*R2)
  D2H (8) = (7.d0*R3*R2 - 15.d0*R4 + 10.d0*R3 - 2.d0*R2)*A
END SUBROUTINE DERIV2_C3_L
                          
SUBROUTINE DERIV_CUBIC (B, A, DH)
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-*
! FIRST DERIVATIVES OF SHAPE FUNCTIONS FOR 1-D
! CUBIC HERMITE ELEMENT   (A C1 ELEMENT)
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-*
Use Precision_Module
 IMPLICIT NONE
 REAL(DP), INTENT(IN)  :: A, B
 REAL(DP), INTENT(OUT) :: DH (4)

! A  = LENGTH OF ELEMENT (SEE SUBR SHAPE_CUBIC)
! B  = COORDINATE OF POINT
! DH = FIRST DERIVATIVES OF H

  DH (1) = 6.d0 * (B * B - B) / A
  DH (2) = 1.d0 - 4.d0 * B + 3.d0 * B * B
  DH (3) = 6.d0 * (B - B * B) / A
  DH (4) = 3.d0 * B * B - 2.d0 * B
END SUBROUTINE DERIV_CUBIC
                       
SUBROUTINE DERIV2_CUBIC (B, A, D2H)
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-*
! SECOND DERIVATIVES OF SHAPE FUNCTIONS FOR 1-D
!  CUBIC HERMITE ELEMENT   (A C1 ELEMENT)
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-*
Use Precision_Module
 IMPLICIT NONE
 REAL(DP), INTENT(IN)  :: A, B
 REAL(DP), INTENT(OUT) :: D2H (4)

! A  = LENGTH OF ELEMENT (SEE SUBR SHAPE_CUBIC)
! B  = COORDINATE OF POINT
! D2H = SECOND DERIVATIVES OF H

  D2H (1) = (12.d0 * B - 6.d0) / A / A
  D2H (2) = (6.d0 * B - 4.d0) / A
  D2H (3) = (6.d0 - 12.d0 * B) / A / A
  D2H (4) = (6.d0 * B - 2.d0) / A
END SUBROUTINE DERIV2_CUBIC
                        
SUBROUTINE DERIV_SHAPE_L_Q_H (NODEDG, LOCATE, NEDGE, N_SPACE, RST,&
                                VALUE, DERIV)
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *
!    SHAPE FUNCTIONS AND DERIVATIVES FOR GENERAL SERENDIPITY
!           LINE, QUAD, OR OR HEXAHEDRON WITH AN
!          ARBITRARY NUMBER OF NODES ON EACH EDGE
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *
Use Precision_Module
 IMPLICIT NONE
 INTEGER, PARAMETER :: MAXDEG = 20
 REAL(DP), INTENT(IN)  :: RST (3)
 INTEGER,  INTENT(IN)  :: NODEDG (12), LOCATE, NEDGE, N_SPACE
 REAL(DP), INTENT(OUT) :: DERIV (3), VALUE

! NODEDG  = NUMBER ON NODES ON 1,4, OR 12 EDGES
! LOCAL   = LOCAL COORDINATE PARALLEL TO EACH EDGE
! NEDGE   = EDGE NUMBER OR CORNER NUMBER OF THE NODE COMPUTED
! N_SPACE = NUMBER OF SPATIAL DIMENSIONS
! RST     = LOCAL COORDINATES FOR EVALUATION
! VALUE   = SHAPE FUNCTION VALUE (RETURNED)
! DERIV   = SHAPE FUNCTION DERIVATIVES (RETURNED)

 REAL(DP) :: BLKCRD (3, 8), POLI2 (3), CDRFN (3), FARSID (3),  &
             CRDEDG (3, MAXDEG + 1)
 INTEGER  :: NEATC (3, 8), NODEOP (2, 12), NODATC (3), LOCAL (12)
 INTEGER  :: ICORD, NSIDE, INODE, NOPV1, NOPV2, ISRFN, &
             INOD1, INOD2, ICOR1, ICOR2
 REAL(DP) :: POLI1, POLI3, CPNUL, PLAN2, DPLA2, DETP2, DETP3, &
             DPOL1, DPOL2, DPOL3
  DATA BLKCRD / -1.d0, -1.d0, -1.d0,    1.d0, -1.d0, -1.d0,   &
                 1.d0,  1.d0, -1.d0,   -1.d0,  1.d0, -1.d0,   &
                -1.d0, -1.d0,  1.d0,    1.d0, -1.d0,  1.d0,   &
                 1.d0,  1.d0,  1.d0,   -1.d0,  1.d0,  1.d0 /
  DATA NEATC /1, 4, 9,   1, 2, 10,   3, 2, 11,   3, 4, 12, &
              5, 8, 9,   5, 6, 10,   7, 6, 11,   7, 8, 12 /
  DATA NODEOP / 7, 8,   8, 5,   5, 6,   6, 7, &
                3, 4,   4, 1,   1, 2,   2, 3, &
                3, 7,   4, 8,   1, 5,   2, 6 /
  DATA LOCAL / -1, -2,  1,  2, -1, -2, &
                1,  2, -3, -3, -3, -3 /
! BLKCRD  = BLOCK CORNER LOCAL COORDINATES
! CRDEDG  = LOCAL COORDINATES OF SIDE NODES JOINING CORNER
! FARSID  = FAR SIDE LOCAL COORDINATE
! LOCATE  = POSITION NUMBER ON EDGE, 0 IF CORNER
! MAXDEG  = MAXIMUM PLOYNOMIAL DEGREE ON ANY SIDE
! NEATC   = THE 1, 2, OR 3 EDGES AT A CORNER
! NODATC  = NUMBER OF SIDE NODES JOINING A CORNER
! NODEOP  = 2 DIAGONALLY OPPOSITE NODES FOR EACH EDGE
!
! VALUE = A(R,S,T)*( P1(R) + P2(S) + P3(T) + CONSTANT )
! DERIV = DA(R,S,T)*( P1(R) + P2(S) + P3(T) + CONSTANT )
!       + A(R,S,T)*( DP1(R) + DP2(S) + DP3(T) )
!
! REF: G. ZAVARISE, ET AL, "AN ALGORITHM FOR GENERATION OF SHAPE
!  FUNCTIONS IN SERENDIPITY ELEMENTS, I.J.N.M.E."
!
!  T:  S      C8 *---E7----* C7   T:  S         8---15----7
!   : /         /.        /:       : /         /.        /:
!   :/         / .       / :       :/         / 22      / :
!   *---R     / E12     /  E11     *---R     /  .      /  20
!           E8   .     E6  :               16   21    /   :
!           /    .    /    :               /    .    /    :
!          /   C4*.../.E3..* C3           /     4.13/.12..3
!         /     .   /     /              /     .   /     /
!     C5 *--E5-----* C6  /              5---------6     /
!        :    .    :    /               :    .    :    11
!        :  E4     :   E2               :  14    19   /
!       E9  .    E10  /                17  .      :  10
!        : .       : /                  : .      18 /
!        :.        :/                   :.        :/
!     C1 *---E1----* C2                 1----9----2
! CORNER NODE & EDGE NUMBERS.   22 NODES: CORNERS, THEN BY EDGES.
!                           CCW IF |T|=1, ELSE IN POSITIVE T.
!                   === 3-D FORM ===
!
!  C4 *---E3----* C3                   4----8----3
!     :    .    :         :S           :         :
!     :         :         :            :         :
!    E4        E2         :            9         :
!     :         :         *---R        :         :
!     :         :                      :         :
!  C1 *---E1----* C2                   1--5-6-7--2
! CORNER NODE & EDGE NUMBERS.  9 NODES: CORNERS, THEN BY EDGE ORDER.
!                    === 2-D FORM ===
!
!  C1 *---E1----* C2                   1--3-4-5--2
! CORNER NODE & EDGE NUMBERS.       5 NODES NUMBERED BY EDGE ORDER.
!                    === 1-D FORM ===
  POLI1 = 1.
  IF (LOCATE == 0) THEN
!
!  SHAPE FUNCTION FOR CORNER NODES
!
    DO ICORD = 1, N_SPACE
      POLI1 = POLI1*(RST (ICORD) + BLKCRD (ICORD, NEDGE))       &
            /(2*BLKCRD (ICORD, NEDGE))
    END DO
    CPNUL = 1.d0
    POLI2 (1) = 0.d0
    POLI2 (2) = 0.d0
    POLI2 (3) = 0.d0
    DO ICORD = 1, N_SPACE
      NSIDE = NEATC (ICORD, NEDGE)
      NODATC (ICORD) = NODEDG (NSIDE) - 2
      IF (NODATC (ICORD) > 0) THEN
        IF (NODATC (ICORD) > MAXDEG) STOP 'ERROR: MAXDEG, DERIV_SHAFN'
        CPNUL = CPNUL - 1.
        POLI2 (ICORD) = 1.
        FARSID (ICORD) = 2.d0/(NODEDG (NSIDE) - 1)
        DO INODE = 1, NODATC (ICORD)
          CRDEDG (ICORD, INODE) = - 1.d0 + FARSID (ICORD)*INODE
          POLI2 (ICORD) = POLI2 (ICORD)*(RST (ICORD)             &
                 - CRDEDG (ICORD, INODE))/(BLKCRD (ICORD, NEDGE) &
                 - CRDEDG (ICORD, INODE))
        END DO
      END IF
    END DO
    VALUE = POLI1*(POLI2 (1) + POLI2 (2) + POLI2 (3) + CPNUL)
  ELSE
!
!  SHAPE FUNCTION FOR EDGE NODES
!
    NOPV1 = NODEOP (1, NEDGE)
    NOPV2 = NODEOP (2, NEDGE)
    ISRFN = ABS (LOCAL (NEDGE))
    FARSID (1) = 2.d0/(NODEDG (NEDGE) - 1)
    CDRFN (1) = - BLKCRD (1, NOPV1)
    CDRFN (2) = - BLKCRD (2, NOPV1)
    CDRFN (3) = - BLKCRD (3, NOPV1)
    CDRFN (ISRFN) = (1.d0 - FARSID (1)*LOCATE)*LOCAL (NEDGE)/ISRFN
    DO ICORD = 1, N_SPACE
      POLI1 = POLI1*(RST (ICORD) - BLKCRD (ICORD, NOPV1))      &
            /(CDRFN (ICORD) - BLKCRD (ICORD, NOPV1))
    END DO
    PLAN2 = (RST (ISRFN) - BLKCRD (ISRFN, NOPV2))              &
          /(CDRFN (ISRFN) - BLKCRD (ISRFN, NOPV2))
    POLI3 = 1.
    NODATC (1) = NODEDG (NEDGE) - 2
    IF (NODATC (1)  > 0) THEN
      IF (NODATC (1)  >  MAXDEG) STOP 'MAXDEG, DERIV_SHAFN'
      DO INODE = 1, NODATC (1)
        CRDEDG (1, INODE) = - 1.d0 + FARSID (1)*INODE
        IF (ABS (CRDEDG (1, INODE) - CDRFN (ISRFN)) > 0.0001) THEN
          POLI3 = POLI3*(RST (ISRFN) - CRDEDG (1, INODE))       &
                /(CDRFN (ISRFN) - CRDEDG (1, INODE))
        END IF
      END DO
    END IF
    VALUE = POLI1*PLAN2*POLI3
  END IF
!
!  DERIVATIVES OF SHAPE FUNCTIONS
!
  DO ICOR1 = 1, N_SPACE
    IF (LOCATE == 0) THEN
!
!  DERIVATIVES FOR CORNER NODES
!
      DPOL1 = POLI2 (1) + POLI2 (2) + POLI2 (3) + CPNUL
      DO ICOR2 = 1, N_SPACE
        IF (ICOR2 /= ICOR1) THEN
          DPOL1 = DPOL1*(RST (ICOR2) + BLKCRD (ICOR2, NEDGE))   &
                /(2*BLKCRD (ICOR2, NEDGE))
        ELSE
          DPOL1 = DPOL1/(2*BLKCRD (ICOR2, NEDGE))
        END IF
      END DO
      DPOL2 = 0.d0
      DO INOD1 = 1, NODATC (ICOR1)
        DETP2 = 1.
        DO INOD2 = 1, NODATC (ICOR1)
          IF (INOD2 /= INOD1) THEN
            DETP2 = DETP2*(RST (ICOR1) - CRDEDG (ICOR1, INOD2)) &
           /(BLKCRD (ICOR1, NEDGE) - CRDEDG (ICOR1, INOD2))
          ELSE
            DETP2 = DETP2/(BLKCRD (ICOR1, NEDGE) - CRDEDG (ICOR1,&
            INOD2))
          END IF
        END DO
        DPOL2 = DPOL2 + DETP2
      END DO
      DPOL2 = DPOL2*POLI1
      DERIV (ICOR1) = DPOL1 + DPOL2
    ELSE
!
!  DERIVATIVES FOR EDGE NODES
!
      DPOL1 = POLI3*PLAN2
      DO ICOR2 = 1, N_SPACE
        IF (ICOR2 /= ICOR1) THEN
          DPOL1 = DPOL1*(RST (ICOR2) - BLKCRD (ICOR2, NOPV1))   &
         /(CDRFN (ICOR2) - BLKCRD (ICOR2, NOPV1))
        ELSE
          DPOL1 = DPOL1/(CDRFN (ICOR2) - BLKCRD (ICOR2, NOPV1))
        END IF
      END DO
      DPLA2 = 0.d0
      DPOL3 = 0.d0
      IF (ICOR1 == ISRFN) THEN
        DPLA2 = POLI1*POLI3/(CDRFN (ISRFN) - BLKCRD (ISRFN,    &
        NOPV2))
        DO INOD1 = 1, NODATC (1)
          IF (ABS (CRDEDG (1, INOD1) - CDRFN (ISRFN))  > 0.0001)  &
          THEN
            DETP3 = 1.
            DO INOD2 = 1, NODATC (1)
              IF (ABS (CRDEDG (1, INOD2) - CDRFN (ISRFN))       &
                  > 0.0001) THEN
                IF (INOD2 /= INOD1) THEN
                  DETP3 = DETP3*(RST (ISRFN) - CRDEDG (1, INOD2))&
                        /(CDRFN (ISRFN) - CRDEDG (1, INOD2))
                ELSE
                  DETP3 = DETP3/(CDRFN (ISRFN) - CRDEDG (1, INOD2))
                END IF
              END IF
            END DO
            DPOL3 = DPOL3 + DETP3
          END IF
        END DO
        DPOL3 = DPOL3*POLI1*PLAN2
      END IF
      DERIV (ICOR1) = DPOL1 + DPLA2 + DPOL3
    END IF
  END DO
END SUBROUTINE DERIV_SHAPE_L_Q_H
                                
SUBROUTINE SHAPE_12_Q (R, S, H)
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-*
! SHAPE FUNCTIONS OF A 12 NODE QUADRILATERAL
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-*
Use Precision_Module
 IMPLICIT NONE
 REAL(DP), INTENT(IN)  :: R, S
 REAL(DP), INTENT(OUT) :: H (12)

 REAL(DP), PARAMETER :: F1 = 1.d0/32, F2 = 9.d0/32
 REAL(DP) :: R_P, R_M, S_P, S_M, R2, S2

! R,S = LOCAL COORDS OF PT     4--11---7---3
! H = ELEM SHAPE FUNCTIONS     I           I
!                              8     S    10
! ELEMENT SKETCH TO RIGHT      I     .R    I
! LOCAL COORD OF NODES:       12           6
! 1-(-1,-1)  3-(+1,+1)         I           I
! SIDES OF ORDER 1, 2, OR 3    1---5---9---2

  R_P = 1 + R ; R_M = 1 - R ; S_P = 1 + S ; S_M = 1 - S
  R2 = R*R ; S2 = S*S
  H ( 1) = F1*R_M*S_M*(-10 + 9*(R2 + S2))
  H ( 2) = F1*R_P*S_M*(-10 + 9*(R2 + S2))
  H ( 3) = F1*R_P*S_P*(-10 + 9*(R2 + S2))
  H ( 4) = F1*R_M*S_P*(-10 + 9*(R2 + S2))
  H ( 5) = F2*(1 - R2)*S_M*(1 - 3*R)
  H ( 6) = F2*(1 - S2)*R_P*(1 - 3*S)
  H ( 7) = F2*(1 - R2)*S_P*(1 + 3*R)
  H ( 8) = F2*(1 - S2)*R_M*(1 + 3*S)
  H ( 9) = F2*(1 - R2)*S_M*(1 + 3*R)
  H (10) = F2*(1 - S2)*R_P*(1 + 3*S)
  H (11) = F2*(1 - R2)*S_P*(1 - 3*R)
  H (12) = F2*(1 - S2)*R_M*(1 - 3*S)
END SUBROUTINE SHAPE_12_Q
                         
SUBROUTINE SHAPE_15_T (R, S, H)
! ******************************************************
! SHAPE FUNCTIONS FOR 15-NODED TRIANGLE
! ******************************************************
Use Precision_Module
 IMPLICIT NONE
 REAL(DP), PARAMETER   :: PT1 = 42.666666666666667d0
 REAL(DP), PARAMETER   :: PT2 = 10.666666666666667d0
 REAL(DP), INTENT(IN)  :: R, S
 REAL(DP), INTENT(OUT) :: H (15)
 REAL(DP) :: CC, T1, T2, T3, T4, T5, T6, T7, T8, T9

!                                   S
!                                   |
!                                   3
! R,S = LOCAL COORDS OF PT          | \
! H = ELEM SHAPE FUNCTIONS          |   \
!                                   6    11
! ELEM_NODES = ELEM INCIDENCES LIST |      \
! ELEMENT SKETCH TO RIGHT           9   15   8
!                                   |          \
! 1@(0,0)  2@(1,0) 2@(0,1)         12   13   14  5
!                                   |             \
!                                   1---4--7---10---2 --R

  CC = 1.d0 - R - S
  T1 = R - 0.25d0   ; T2 = R - 0.5d0   ; T3 = R - 0.75d0
  T4 = S - 0.25d0   ; T5 = S - 0.5d0   ; T6 = S - 0.75d0
  T7 = CC - 0.25d0  ; T8 = CC - 0.5d0  ; T9 = CC - 0.75d0
  H (1) = PT2*CC*T7*T8*T9
  H (2) = PT2*R*T1*T2*T3
  H (3) = PT2*S*T4*T5*T6
  H (4) = PT1*CC*R*T7*T8
  H (5) = PT1*R*S*T1*T2
  H (6) = PT1*S*CC*T4*T5
  H (7) = 64.d0*CC*R*T1*T7
  H (8) = 64.d0*R*S*T1*T4
  H (9) = 64.d0*S*CC*T4*T7
  H (10) = PT1*CC*R*T1*T2
  H (11) = PT1*R*S*T4*T5
  H (12) = PT1*S*CC*T7*T8
  H (13) = 128.d0*R*S*CC*T7
  H (14) = 128.d0*R*S*T1*CC
  H (15) = 128.d0*R*S*CC*T4
END SUBROUTINE SHAPE_15_T
                         
SUBROUTINE SHAPE_16_Q (R, S, H)
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-*
! LAGRANGIAN 16 NODE QUAD IN NATURAL COORDINATES
!   USING TENSOR PRODUCTS OF 1D BASIS
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-*
 Use Precision_Module
 IMPLICIT NONE
 REAL(DP), INTENT(IN)  :: R, S
 REAL(DP), INTENT(OUT) :: H (16)
 REAL(DP)              :: HR (4), HS (4)

! 4   4--11---7---3   The active Lagrangian Q16 elements
! .   I     S     I   have nodes input in the order shown
! 3   8  16 . 15 10   to the left.
! .   I     ..R   I
! 2  12  13   14  6
! .   I           I
! 1   1---5---9---2
!
!     1---2---3---4 -->R

! Evaluate the 1D interpolations
  CALL SHAPE_4_L (R, HR)
  CALL SHAPE_4_L (S, HS)

! Form tensor products
  H (1) = HR (1)*HS (1)
  H (2) = HR (4)*HS (1)
  H (3) = HR (4)*HS (4)
  H (4) = HR (1)*HS (4)
  H (5) = HR (2)*HS (1)
  H (6) = HR (4)*HS (2)
  H (7) = HR (3)*HS (4)
  H (8) = HR (1)*HS (3)
  H (9) = HR (3)*HS (1)
  H (10) = HR (4)*HS (3)
  H (11) = HR (2)*HS (4)
  H (12) = HR (1)*HS (2)
  H (13) = HR (2)*HS (2)
  H (14) = HR (3)*HS (2)
  H (15) = HR (3)*HS (3)
  H (16) = HR (2)*HS (3)
END SUBROUTINE SHAPE_16_Q

SUBROUTINE SHAPE_16_QS (R, S, H)
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-**
!   SHAPE FUNCTIONS FOR SERENDIPITY QUAD WITH 16 NODES
!                A BI-4TH ORDER ELEMENT
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-**
Use Precision_Module
 IMPLICIT NONE
 REAL(DP), PARAMETER   :: PT667 = 0.66666666666667d0
 REAL(DP), INTENT(IN)  :: R, S
 REAL(DP), INTENT(OUT) :: H (16)
 REAL(DP) :: RR, SS, RS, RP, RM, S_P, SM

! R,S = LOCAL COORDS OF PT          4--15--11---7---3
! H = ELEM SHAPE FUNCTIONS          I               I
!                                   8       S      14
!                                   I       .       I
! ELEMENT SKETCH TO RIGHT          12       +..R   10
!                                   I               I
! 1@(-1,-1)  3@(+1,+1)             16               6
!                                   I               I
!                                   1---5---9--13---2

  RR = R*R     ; SS = S*S    ; RS = R*S
  RP = 1.d0 + R  ; RM = 1.d0 - R
  S_P = 1.d0 + S  ; SM = 1.d0 - S
  H (1) = RM*SM*(-R*(4.d0*RR - 1.d0) - S*(4.d0*SS - 1.d0) - 3.d0)/12.d0
  H (5) = - PT667*R*SM*RM*RP*(1.d0 - 2.d0*R)
  H (9) = 0.5d0*RM*RP*(1.d0 - 4.d0*RR)*SM
  H (13) = PT667*R*SM*RM*RP*(1.d0 + 2.d0*R)
  H (2) = RP*SM*(R*(4.d0*RR - 1.d0) - S*(4.d0*SS - 1.d0) - 3.d0)/12.d0
  H (6) = - PT667*S*RP*SM*S_P*(1.d0 - 2.d0*S)
  H (10) = 0.5d0*SM*S_P*(1.d0 - 4.d0*SS)*RP
  H (14) = PT667*S*RP*SM*S_P*(1.d0 + 2.d0*S)
  H (3) = RP*S_P*(R*(4.d0*RR - 1.d0) + S*(4.d0*SS - 1.d0) - 3.d0)/12.d0
  H (7) = PT667*R*S_P*RM*RP*(1.d0 + 2.d0*R)
  H (11) = 0.5d0*RM*RP*(1.d0 - 4.d0*RR)*S_P
  H (15) = - PT667*R*S_P*RM*RP*(1.d0 - 2.d0*R)
  H (4) = RM*S_P*(-R*(4.d0*RR - 1.d0) + S*(4.d0*SS - 1.d0) - 3.d0)/12.d0
  H (8) = PT667*S*SM*S_P*(1.d0 + 2.d0*S)*RM
  H (12) = 0.5d0*SM*S_P*(1.d0 - 4.d0*SS)*RM
  H (16) = - PT667*S*RM*SM*S_P*(1.d0 - 2.d0*S)
END SUBROUTINE SHAPE_16_QS
                          
SUBROUTINE SHAPE_16_R (R, S, A, B, H)
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-*
! WARNING, this code is not fully checked. No known bugs.
! C1 RECTANGULAR ELEMENT IN UNIT COORDINATES
!   USING TENSOR PRODUCTS OF 1D BASIS
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-*
Use Precision_Module
 IMPLICIT NONE
 REAL(DP), INTENT(IN)  :: R, S, A, B
 REAL(DP), INTENT(OUT) :: H (16)
 REAL(DP)              :: HR (4), HS (4)

! DOF ARE W W,X W,Y W,XY AT EACH NODE (N_G_DOF=4)
! X // R, Y // S.              S
! A = PHYSICAL LENGTH IN X     4 -------- 3
! B = PHYSICAL LENGTH IN Y     I          I
! R,S = LOCAL UNIT COORDS      I          I
! 1@(0,0), 3@(1,1)             1 -------- 2 ->R

! Evaluate the 1D interpolations
  CALL SHAPE_C1_L (R, A, HR)
  CALL SHAPE_C1_L (S, B, HS)

! Form tensor products
  H (1) = HR (1) * HS (1)
  H (2) = HR (2) * HS (1)
  H (3) = HR (1) * HS (2)
  H (4) = HR (2) * HS (2)
  H (5) = HR (3) * HS (1)
  H (6) = HR (4) * HS (1)
  H (7) = HR (3) * HS (2)
  H (8) = HR (4) * HS (2)
  H (9) = HR (3) * HS (3)
  H (10) = HR (4) * HS (3)
  H (11) = HR (3) * HS (4)
  H (12) = HR (4) * HS (4)
  H (13) = HR (1) * HS (3)
  H (14) = HR (2) * HS (3)
  H (15) = HR (1) * HS (4)
  H (16) = HR (2) * HS (4)
END SUBROUTINE SHAPE_16_R
                         
SUBROUTINE SHAPE_17_Q (R, S, H)
! ******************************************************************
!   SHAPE FUNCTIONS FOR A SERENDIPITY QUAD WITH 17 NODES
! ******************************************************************
Use Precision_Module
 IMPLICIT NONE
 REAL(DP), PARAMETER   :: PT667 = 0.6666666666666667d0
 REAL(DP), INTENT(IN)  :: R, S
 REAL(DP), INTENT(OUT) :: H (17)
 REAL(DP)              :: RR, SS, RS, RP, RM, S_P, SM

! R,S = LOCAL COORDS OF PT          4--15--11---7---3
! H = ELEM SHAPE FUNCTIONS          I               I
!                                   8       S      14
!                                   I       I       I
! ELEMENT SKETCH TO RIGHT          12      17-R    10
!                                   I               I
! 1@(-1,-1)  3@(+1,+1)             16               6
! 17@(0,0)                          I               I
!                                   1---5---9--13---2

  RR = R*R  ; SS = S*S  ; RS = R*S
  RP = 1.d0 + R ; RM = 1.d0 - R
  S_P = 1.d0 + S ; SM = 1.d0 - S
  H (1) = RM*SM*( - 4.d0*R*(RR - 1.d0) - 4.d0*S*(SS - 1.d0) + 3.d0*RS)/12.d0
  H (5) = - PT667*R*SM*RM*RP*(1.d0 - 2.d0*R)
  H (9) = 0.5d0*RM*RP*( - S - 4.d0*RR)*SM
  H (13) = PT667*R*SM*RM*RP*(1.d0 + 2.d0*R)
  H (2) = RP*SM*(4.d0*R*(RR - 1.d0) - 4.d0*S*(SS - 1.d0) - 3.d0*RS)/12.d0
  H (6) = - PT667*S*RP*SM*S_P*(1.d0 - 2.d0*S)
  H (10) = 0.5d0*SM*S_P*(R - 4.d0*SS)*RP
  H (14) = PT667*S*RP*SM*S_P*(1.d0 + 2.d0*S)
  H (3) = RP*S_P*(4.d0*R*(RR - 1.d0) + 4.d0*S*(SS - 1.d0) + 3.d0*RS)/12.d0
  H (7) = PT667*R*S_P*RM*RP*(1.d0 + 2.d0*R)
  H (11) = 0.5d0*RM*RP*(S - 4.d0*RR)*S_P
  H (15) = - PT667*R*S_P*RM*RP*(1.d0 - 2.d0*R)
  H (4) = RM*S_P*( - 4.d0*R*(RR - 1.d0) + 4.d0*S*(SS - 1.d0) - 3.d0*RS)/12.d0
  H (8) = PT667*S*SM*S_P*(1.d0 + 2.d0*S)*RM
  H (12) = 0.5d0*SM*S_P*( - R - 4.d0*SS)*RM
  H (16) = - PT667*S*RM*SM*S_P*(1.d0 - 2.d0*S)
  H (17) = RM*RP*SM*S_P
END SUBROUTINE SHAPE_17_Q
                         
SUBROUTINE SHAPE_8_20_H (R, S, T, H, ELEM_NODES)
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-*
! ELEMENT INTERPOLATION FUNCTIONS FOR AN 8 TO 20 NODE HEXAHEDRON
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-*
Use Precision_Module
 IMPLICIT NONE
 REAL(DP), INTENT(IN)  :: R, S, T
 REAL(DP), INTENT(OUT) :: H (20)
 INTEGER, INTENT(IN)  :: ELEM_NODES (20)
 INTEGER              :: I1 (20), I2 (20), K, K1, K2
 REAL(DP)                 :: RP, RM, S_P, SM, TP, TM, RZ, SZ, TZ, HK
  DATA I1/8*0, 1, 2, 3, 4, 5, 6, 7, 8, 1, 2, 3, 4/
  DATA I2/8*0, 2, 3, 4, 1, 6, 7, 8, 5, 5, 6, 7, 8/

! R,S,T = LOCAL COORDINATES OF THE POINT   -1 LE (R,S,T) LE +1
! H     = ELEMENT INTERPOLATION FUNCTIONS
!         H(I) = 0 IF ELEM_NODES(I) = 0
! ELEM_NODES = ARRAY OF ELEMENT INCIDENCES,
!  IF ELEM_NODES(I)=0 THEN LOCAL NODE I IS NOT CONSIDERED IN ANALYSIS
! I1, I2 = CORNER NODES OF TWELVE EDGES
!
! A SKETCH OF THE LOCAL NODES...        T          3 *----*----* 2
!                                       :          / .    10  /:
! FACES DEFINED BY R.H.R. ABOUT         :         /  .       / :
!    THE POSITIVE LOCAL AXES            *---S    /   *19    / *18
!  +R,(8,5,1,4,16,17,12,20)           /       11*    .    9*   :
!  -R,(7,6,2,3,14,18,10,19)          /         /     .    /    :
!  +S,(6,2,1,5,18,9,17,13)           R        /     7*.../*....* 6
!  -S,(7,3,4,8,19,11,20,15)                  /   12 .   /14   /
!  +T,(3,4,1,2,11,12,9,10)                 4 *----*----*1    /
!  -T,(7,8,5,6,15,16,13,14)                  :    .    :    /
! AND  R             FACE                    :   *15   :   *13
!     7 V           7    V                20 *  .    17*  /
!    7   V         7      V                  : .       : /
!   T <<< S   COORD2 <<<< COORD1             :.   16   :/
! ARE THE FACE COORD. PERMUTATIONS         8 *----*----* 5

  RP = 0.5d0*(1.d0 + R) ; S_P = 0.5d0*(1.d0 + S) ; TP = 0.5d0*(1.d0 + T)
  RM = 0.5d0*(1.d0 - R) ; SM = 0.5d0*(1.d0 - S) ; TM = 0.5d0*(1.d0 - T)
  RZ = 1.d0 - R*R     ; SZ = 1.d0 - S*S     ; TZ = 1.d0 - T*T
  H (1) = TP*S_P*RP
  H (2) = TP*S_P*RM
  H (3) = TP*SM*RM
  H (4) = TP*SM*RP
  H (5) = TM*S_P*RP
  H (6) = TM*S_P*RM
  H (7) = TM*SM*RM
  H (8) = TM*SM*RP

!  QUADRATIC EDGE BUBBLES
  H (9) = TP*S_P*RZ*0.5d0
  H (10) = TP*SZ*RM*0.5d0
  H (11) = TP*SM*RZ*0.5d0
  H (12) = TP*SZ*RP*0.5d0
  H (13) = TM*S_P*RZ*0.5d0
  H (14) = TM*SZ*RM*0.5d0
  H (15) = TM*SM*RZ*0.5d0
  H (16) = TM*SZ*RP*0.5d0
  H (17) = TZ*S_P*RP*0.5d0
  H (18) = TZ*S_P*RM*0.5d0
  H (19) = TZ*SM*RM*0.5d0
  H (20) = TZ*SM*RP*0.5d0

!  LOOP OVER TWELVE ELEMENT EDGES
  DO K = 9, 20
    IF ( ELEM_NODES (K)  == 0 ) THEN
!      SET UNUSED EDGE BUBBLE TO ZERO
      H (K) = 0.d0
    ELSE
!      ENRICH THE TWO CORNERS ON THE EDGE
      HK = H (K)
      K1 = I1 (K)
      K2 = I2 (K)
      H (K1) = H (K1) - HK
      H (K2) = H (K2) - HK
      H (K) = HK + HK
    END IF
  END DO
END SUBROUTINE SHAPE_8_20_H
                           
SUBROUTINE SHAPE_20_H (R, S, T, H)
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-*
!   ELEMENT INTERPOLATION FUNCTIONS FOR A 20 NODE HEXAHEDRON
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-*
Use Precision_Module
 IMPLICIT NONE
 REAL(DP), INTENT(IN)  :: R, S, T
 REAL(DP), INTENT(OUT) :: H (20)
 INTEGER              :: I1 (20), I2 (20), K, K1, K2
 REAL(DP)                 :: RP, RM, S_P, SM, TP, TM, RZ, SZ, TZ, HK
  DATA I1/8*0, 1, 2, 3, 4, 5, 6, 7, 8, 1, 2, 3, 4/
  DATA I2/8*0, 2, 3, 4, 1, 6, 7, 8, 5, 5, 6, 7, 8/

! R,S,T  = LOCAL COORDINATES OF THE POINT   -1 LE (R,S,T) LE +1
! H      = ELEMENT INTERPOLATION FUNCTIONS
! I1, I2 = CORNER NODES OF TWELVE EDGES
!
! A SKETCH OF THE LOCAL NODES...        T          3 *----*----* 2
!                                       :          / .    10  /:
! FACES DEFINED BY R.H.R. ABOUT         :         /  .       / :
!    THE POSITIVE LOCAL AXES            *---S    /   *19    / *18
!  +R,(8,5,1,4,16,17,12,20)           /       11*    .    9*   :
!  -R,(7,6,2,3,14,18,10,19)          /         /     .    /    :
!  +S,(6,2,1,5,18,9,17,13)           R        /     7*.../*....* 6
!  -S,(7,3,4,8,19,11,20,15)                  /   12 .   /14   /
!  +T,(3,4,1,2,11,12,9,10)                 4 *----*----*1    /
!  -T,(7,8,5,6,15,16,13,14)                  :    .    :    /
! AND  R             FACE                    :   *15   :   *13
!     7 V           7    V                20 *  .    17*  /
!    7   V         7      V                  : .       : /
!   T <<< S   COORD2 <<<< COORD1             :.   16   :/
! ARE THE FACE COORD. PERMUTATIONS         8 *----*----* 5

  RP = 0.5d0*(1.d0 + R) ; S_P = 0.5d0*(1.d0 + S) ; TP = 0.5d0*(1.d0 + T)
  RM = 0.5d0*(1.d0 - R) ; SM = 0.5d0*(1.d0 - S) ; TM = 0.5d0*(1.d0 - T)
  RZ = 1.d0 - R*R     ; SZ = 1.d0 - S*S     ; TZ = 1.d0 - T*T
  H (1) = TP*S_P*RP
  H (2) = TP*S_P*RM
  H (3) = TP*SM*RM
  H (4) = TP*SM*RP
  H (5) = TM*S_P*RP
  H (6) = TM*S_P*RM
  H (7) = TM*SM*RM
  H (8) = TM*SM*RP

!  QUADRATIC EDGE BUBBLES
  H (9) = TP*S_P*RZ*0.5d0
  H (10) = TP*SZ*RM*0.5d0
  H (11) = TP*SM*RZ*0.5d0
  H (12) = TP*SZ*RP*0.5d0
  H (13) = TM*S_P*RZ*0.5d0
  H (14) = TM*SZ*RM*0.5d0
  H (15) = TM*SM*RZ*0.5d0
  H (16) = TM*SZ*RP*0.5d0
  H (17) = TZ*S_P*RP*0.5d0
  H (18) = TZ*S_P*RM*0.5d0
  H (19) = TZ*SM*RM*0.5d0
  H (20) = TZ*SM*RP*0.5d0

!  LOOP OVER TWELVE ELEMENT EDGES
  DO K = 9, 20
!    ENRICH THE TWO CORNERS ON THE EDGE
    HK = H (K)
    K1 = I1 (K)
    K2 = I2 (K)
    H (K1) = H (K1) - HK
    H (K2) = H (K2) - HK
    H (K) = HK + HK
  END DO
END SUBROUTINE SHAPE_20_H
                           
SUBROUTINE SHAPE_2_L (R, H)
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-**
! SHAPE FUNCTIONS OF A 2 NODE LINE ELEMENT
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-**
Use Precision_Module
 IMPLICIT NONE
 REAL(DP), INTENT(IN)  :: R
 REAL(DP), INTENT(OUT) :: H (2)

! R IS NATURAL COORD.   R=-1  1------------2  R=1

  H (1) = 0.5d0*(1.d0 - R)
  H (2) = 0.5d0*(1.d0 + R)
END SUBROUTINE SHAPE_2_L
                        
SUBROUTINE SHAPE_3_T (S, T, H)
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-*
! SHAPE FUNCTIONS FOR A THREE NODE UNIT TRIANGLE
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-*
Use Precision_Module
 IMPLICIT NONE
 REAL(DP), INTENT(IN)  :: S, T
 REAL(DP), INTENT(OUT) :: H (3)

! S,T = LOCAL COORDINATES OF THE POINT    3     T
! H   = SHAPE FUNCTIONS                   . .   .
! NODAL COORDS 1-(0,0)  2-(1,0)  3-(0,1)  1..2  0..S

  H (1) = 1.d0 - S - T
  H (2) = S
  H (3) = T
END SUBROUTINE SHAPE_3_T
                        
!SUBROUTINE SHAPE_4_12_Q (R, S, H, ELEM_NODES)
!! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-*
!!    SHAPE FUNCTIONS OF 4 TO 12 NODE QUADRILATERAL
!! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-*
!Use Precision_Module
! IMPLICIT NONE
! REAL(DP), INTENT(IN)  :: R, S
! REAL(DP), INTENT(OUT) :: H (12)
! INTEGER, INTENT(IN)  :: ELEM_NODES (12)
! INTEGER              :: IP (4), JP (4), NEXT (4), I, K
! REAL(DP)                 :: P, Q, TEMP
!  DATA IP, JP, NEXT/1, 0, - 1, 0, 0, 1, 0, - 1, 2, 3, 4, 1/
!
!! R,S = LOCAL COORDS OF PT          4--11---7---3
!! H = ELEM SHAPE FUNCTIONS          :     S     :
!! ELEM_NODES = ELEM INCIDENCES LIST 8     :    10
!! ELEMENT SKETCH TO RIGHT           :     *..R  :
!! LOCAL COORD OF NODES:            12           6
!! 1-(-1,-1)  3-(+1,+1)              :           :
!! SIDES OF ORDER 1, 2, OR 3         1---5---9---2
!
!!-->   GENERATE FOUR NODE BILINEAR QUADRILATERAL
!  CALL SHAPE_4_Q (R, S, H)
!
!!-->   LOOP OVER SIDES
!  DO I = 1, 4
!    H (I + 4) = 0.d0
!    H (I + 8) = 0.d0
!!    IS SIDE HIGHER THAN LINEAR ?
!    IF ( ELEM_NODES (I + 4)  > 0 ) THEN
!      K = NEXT (I)
!!      FIND PT RELATIVE TO SIDE (CYCLIC COORDINATES)
!      P = R*IP (I) + S*JP (I)
!      Q = - R*JP (I) + S*IP (I)
!      TEMP = (1.d0 - Q)*0.5d0
!!-->       IS SIDE QUADRATIC OR CUBIC ?
!      IF ( ELEM_NODES (I + 8)  > 0 ) THEN
!!        CUBIC
!        H (I + 8) = TEMP*(1.d0 - 3.d0*P - P*P + 3.d0*P*P*P)*9.d0/16.d0
!        H (I + 4) = TEMP*(1.d0 + 3.d0*P - P*P - 3.d0*P*P*P)*9.d0/16.d0
!!        CORRECT CORNER POINTS FOR CUBIC
!        H (I) = H (I) - H (I + 4)*2.d0/3.d0 - H (I + 8)/3.d0
!        H (K) = H (K) - H (I + 8)*2.d0/3.d0 - H (I + 4)/3.d0
!      ELSE
!!        QUADRATIC
!        H (I + 4) = TEMP*(1.d0 - P*P)
!!        CORRECT CORNER POINTS FOR QUADRATIC
!        H (I) = H (I) - H (I + 4)*0.5d0
!        H (K) = H (K) - H (I + 4)*0.5d0
!      END IF
!    END IF
!  END DO
!END SUBROUTINE SHAPE_4_12_Q
                           
SUBROUTINE SHAPE_4_L (X, H)
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-*
! CALCULATE SHAPE FUNCTION OF A 4 NODE LINE ELEMENT
!         IN NATURAL COORDINATES
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-*
Use Precision_Module
 IMPLICIT NONE
 REAL(DP), INTENT(IN)  :: X
 REAL(DP), INTENT(OUT) :: H (4)

! H = ELEMENT SHAPE FUNCITON
! X = LOCAL COORDIATE OF POINT, -1 <= X <= +1
! LOCAL NODE COORD. ARE    1----2---3----4

  H (1) = (1.d0 - X)*(3.d0*X + 1.d0)*(3.d0*X - 1.d0)/16.d0
  H (2) = 9.d0*(1.d0 + X)*(X - 1.d0)*(3.d0*X - 1.d0)/16.d0
  H (3) = 9.d0*(1.d0 + X)*(1.d0 - X)*(3.d0*X + 1.d0)/16.d0
  H (4) = (1.d0 + X)*(3.d0*X + 1.d0)*(3.d0*X - 1.d0)/16.d0
END SUBROUTINE SHAPE_4_L
                        
SUBROUTINE SHAPE_4_P (R, S, T, H)
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-*
! SHAPE FUNCTIONS OF A 4 NODE TETRAHEDRON ELEMENT
!          (PYRAMID) IN UNIT COORDINATES
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-*
Use Precision_Module               !     s
 IMPLICIT NONE                     !     |
 REAL(DP), INTENT(IN)  :: R, S, T  !     3
 REAL(DP), INTENT(OUT) :: H (4)    !     | .
                                   !   . 1---2 --r
  H (1) = 1.d0 - R - S - T         !    /  .
  H (2) = R                        !   4
  H (3) = S                        !  /
  H (4) = T                        ! t
END SUBROUTINE SHAPE_4_P
                        
SUBROUTINE SHAPE_4_Q (R, S, H)
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-*
! SHAPE FUNCTIONS OF A 4 NODE PARAMETRIC QUAD
!          IN NATURAL COORDINATES
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-*
Use Precision_Module
 IMPLICIT NONE
 REAL(DP), INTENT(IN)  :: R, S
 REAL(DP), INTENT(OUT) :: H (4)
 REAL(DP)              :: RP, RM, S_P, SM

! (R,S) = A POINT IN THE NATURAL COORDS     4--3
! H     = LOCAL INTERPOLATION FUNCTIONS     I  I
! H(I)  = 0.25d0*(1+R*R(I))*(1+S*S(I))      I  I
! R(I)  = LOCAL R-COORDINATE OF NODE I      1--2
! LOCAL COORDS, 1=(-1,-1)   3=(+1,+1)

  RP = 1.d0 + R ; RM = 1.d0 - R
  S_P = 1.d0 + S ; SM = 1.d0 - S
  H (1) = 0.25d0*RM*SM
  H (2) = 0.25d0*RP*SM
  H (3) = 0.25d0*RP*S_P
  H (4) = 0.25d0*RM*S_P
END SUBROUTINE SHAPE_4_Q
                        
SUBROUTINE SHAPE_4_R (R, S, H)
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-*
! C0 RECTANGULAR ELEMENT IN UNIT COORDINATES
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-*
Use Precision_Module
 IMPLICIT NONE
 REAL(DP), INTENT(IN)  :: R, S
 REAL(DP), INTENT(OUT) :: H (4)

! DOF ARE W AT EACH NODE (N_G_DOF=1)
! X // R, Y // S.              S
! A = PHYSICAL LENGTH IN X     4 -------- 3
! B = PHYSICAL LENGTH IN Y     I          I
! R,S = LOCAL UNIT COORDS      I          I
! 1@(0,0), 3@(1,1)             1 -------- 2 ->R

  H (1) = 1.d0 - R - S + R*S
  H (2) = R - R*S
  H (3) = R*S
  H (4) = S - R*S
END SUBROUTINE SHAPE_4_R
                        
SUBROUTINE SHAPE_4_T (R, S, H)
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-**
! UNIT TRIANGLE WITH CENTROID (4TH) NODE
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-**
Use Precision_Module
 IMPLICIT NONE
 REAL(DP), INTENT(IN)  :: R, S
 REAL(DP), INTENT(OUT) :: H (4)

! R,S = UNIT COORDINATES, SEE SHAPE_4_T
! H = INTERPOLATION FUNCTIONS
! NODE 1@(0,0), 2@(1,0), 3@(0,1), 4@(1/3,1/3)

  H (4) = 27.d0*R*S*(1.d0 - R - S)
  H (1) = 1.d0 - R - S - H (4)/3.d0
  H (2) = R - H (4)/3.d0
  H (3) = S - H (4)/3.d0
END SUBROUTINE SHAPE_4_T
                        
SUBROUTINE SHAPE_7_T (R, S, H)
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-**
! QUADRATIC UNIT TRIANGLE PLUS CENTROID NODE
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-**
Use Precision_Module
 IMPLICIT NONE
 REAL(DP), PARAMETER   :: PT44 = 0.44444444444444d0
 REAL(DP), INTENT(IN)  :: R, S
 REAL(DP), INTENT(OUT) :: H (7)

! R,S = UNIT COORDINATES, SEE SHAPE_6_T
! H = INTERPOLATION FUNCTIONS
! NODE 7@(1/3,1/3), OTHERS LIKE SHAPE_6_T

  H (7) = 27.d0*R*S*(1.d0 - R - S)
  H (1) = 1.d0 - 3.d0*R - 3.d0*S + 2.d0*R*R + 4.d0*R*S + 2.d0*S*S + H (7)/9.d0
  H (2) = 2.d0*R*R - R + H (7)/9.d0
  H (3) = 2.d0*S*S - S + H (7)/9.d0
  H (4) = 4.d0*(R - R*R - R*S) - H (7)*PT44
  H (5) = 4.d0*R*S - H (7)*PT44
  H (6) = 4.d0*(S - R*S - S*S) - H (7)*PT44
END SUBROUTINE SHAPE_7_T
                        
SUBROUTINE SHAPE_8_H (R, S, T, H)
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-*
! SHAPE FUNCTIONS OF 8 NODE PARAMETRIC HEXAHEDRON
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-*
Use Precision_Module
 IMPLICIT NONE
 REAL(DP), INTENT(IN)  :: R, S, T
 REAL(DP), INTENT(OUT) :: H (8)
 REAL(DP)              :: RP, RM, S_P, SM, TP, TM

! R,S,T = LOCAL COORDS OF PT !     R         3 *---------* 2
! H = ELEM SHAPE FUNCTIONS   !     |         / .        /|
! NODES ORDERED BY RHR       !     |        /  .       / |
! ABOUT THE R-AXIS           !     *---T   /   .      /  |
! LOCAL COORD:1=(1,1,1)      !   /        /    .     /   |
! 4=(1,1,-1)  7=(-1,-1,-1)   !  /        /     .    /    |
                             !  S       /     7*.../.....* 6
 RP = 1.d0 + R               !         /      .   /     /
 RM = 1.d0 - R               !       4 *---------*1    /
 S_P = 1.d0 + S              !         |    .    |    /
 SM = 1.d0 - S               !         |   .     |   /
 TP = 1.d0 + T               !         |  .      |  /
 TM = 1.d0 - T               !         | .       | /
                             !         |.        |/
                             !       8 *---------* 5
  H (1) = 0.125d0*RP*S_P*TP
  H (2) = 0.125d0*RP*SM*TP
  H (3) = 0.125d0*RP*SM*TM
  H (4) = 0.125d0*RP*S_P*TM
  H (5) = 0.125d0*RM*S_P*TP
  H (6) = 0.125d0*RM*SM*TP
  H (7) = 0.125d0*RM*SM*TM
  H (8) = 0.125d0*RM*S_P*TM
END SUBROUTINE SHAPE_8_H
                        
SUBROUTINE SHAPE_9_Q (R, S, H)
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-*
! SHAPE FUNCTIONS FOR 9-NODED QUAD
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-*
Use Precision_Module
 IMPLICIT NONE
 REAL(DP), INTENT(IN)  :: R, S
 REAL(DP), INTENT(OUT) :: H (9)
 REAL(DP)              :: RP, RM, S_P, SM

! R,S = LOCAL COORDS OF PT          4-----7-----3
! H = ELEM SHAPE FUNCTIONS          I     S     I
!                                   I     .     I
! ELEMENT SKETCH TO RIGHT           8     9..R  6
! 1@(-1,-1)  3@(+1,+1)              I           I
!                                   I           I
!                                   1-----5-----2

  RM = R - 1.D0 ; SM = S - 1.D0
  RP = R + 1.D0 ; S_P = S + 1.D0
  H (1) = 0.25D0*S*SM*R*RM
  H (2) = 0.25D0*S*SM*R*RP
  H (3) = 0.25D0*S*S_P*R*RP
  H (4) = 0.25D0*S*S_P*R*RM
  H (5) = - 0.5D0*S*SM*RP*RM
  H (6) = - 0.5D0*S_P*SM*R*RP
  H (7) = - 0.5D0*S*S_P*RP*RM
  H (8) = - 0.5D0*S_P*SM*R*RM
  H (9) = S_P*SM*RP*RM
END SUBROUTINE SHAPE_9_Q
                        
SUBROUTINE SHAPE_C0_L (R, H)
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-**
! C0 LINE ELEM IN UNIT COORDINATES
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-**
Use Precision_Module
 IMPLICIT NONE
 REAL(DP), INTENT(IN)  :: R
 REAL(DP), INTENT(OUT) :: H (2)

! H = SHAPE FUNCTIONS (R=0) 1----2 (R=1) ->R
! DOF ARE FUNCTION VALUES AT EACH NODE

  H (1) = 1.d0 - R
  H (2) = R
END SUBROUTINE SHAPE_C0_L
                         
SUBROUTINE SHAPE_C1_L (R, A, H)
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-*
! SHAPE FUNCTIONS FOR CUBIC HERMITE IN UNIT COORDINATES
!                    ( A C1 ELEMENT )
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-*
Use Precision_Module
 IMPLICIT NONE
 REAL(DP), INTENT(IN)  :: R, A
 REAL(DP), INTENT(OUT) :: H (4)

! A = PHYSICAL LENGTH OF ELEMENT     1----------2 ---> R
! R = LOCAL COORDINATE OF POINT     R=0         R=1
! H = SHAPE FUNCTIONS ARRAY
! DOF ARE FUNCTION AND SLOPE, WRT X, AT EACH NODE (N_G_DOF=2)
! D()/DX = D()/DR DR/DX = 1/A * D()/DR

  H(1) = 1.d0 - 3.0*R*R + 2.0*R*R*R
  H(2) = (R - 2.0*R*R + R*R*R)*A
  H(3) =  3.0*R*R - 2.0*R*R*R
  H(4) = (R*R*R - R*R)*A
END SUBROUTINE SHAPE_C1_L
                         
SUBROUTINE SHAPE_C1_L2_NC (R, A, H)
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-*
! SHAPE FUNCTIONS FOR CUBIC HERMITE IN NATURAL COORDINATES
!                    ( A C1 ELEMENT )
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-*
Use Precision_Module
 IMPLICIT NONE
 REAL(DP), INTENT(IN)  :: R, A
 REAL(DP), INTENT(OUT) :: H (4)

! A = PHYSICAL LENGTH OF ELEMENT     1-----+-----2 ---> R
! R = LOCAL COORDINATE OF POINT     R=-1   0    R=1
! H = SHAPE FUNCTIONS ARRAY
! DOF ARE FUNCTION AND SLOPE, WRT X, AT EACH NODE (N_G_DOF=2)
! D()/DX = D()/DR DR/DX = 2/A * D()/DR

  H(1) = (2 - 3*R + R*R*R)/4        
  H(2) = (1 - R - R*R + R*R*R)*A/8
  H(3) = (2 + 3*R - R*R*R)/4    
  H(4) = (-1 - R + R*R + R*R*R)*A/8
END SUBROUTINE SHAPE_C1_L2_NC
                         
SUBROUTINE SHAPE_C2_L (R, A, H)
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-**
! SHAPE FUNCTIONS FOR FIFTH ORDER HERMITE LINE ELEMENT
!            ( A C2 ELEMENT IN UNIT COORDINATES)
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-**
Use Precision_Module
 IMPLICIT NONE
 REAL(DP), INTENT(IN)  :: R, A
 REAL(DP), INTENT(OUT) :: H (6)
 REAL(DP)              :: R3

! A = PHYSICAL LENGTH OF ELEMENT     1----------2 -> R
! R = LOCAL COORDINATE OF POINT     R=0         R=1
! H = SHAPE FUNCTIONS ARRAY
! DOF ARE VALUE, SLOPE, CURVATURE AT EACH END (WRT X)
! D()/DX = D()/DR DR/DX = 1/A*D()/DR

  R3 = R*R*R
  H (1) = 1.d0 - 10.d0*R3 + 15.d0*R3*R - 6.d0*R3*R*R
  H (2) = (R - 6.d0*R3 + 8.d0*R3*R - 3.d0*R3*R*R)*A
  H (3) = (R*R - 3.d0*R3 + 3.d0*R3*R - R3*R*R)*A*A/2.d0
  H (4) = 10.d0*R3 - 15.d0*R3*R + 6.d0*R3*R*R
  H (5) = (7.d0*R3*R - 3.d0*R3*R*R - 4.d0*R3)*A
  H (6) = (R3 - 2.d0*R3*R + R3*R*R)*A*A/2.d0
END SUBROUTINE SHAPE_C2_L
                         
SUBROUTINE SHAPE_C3_L (R, A, H)
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-*
! SHAPE FUNCTIONS FOR SEVENTH ORDER HERMITE LINE ELEMENT
!            ( A C3 ELEMENT IN UNIT COORDINATES)
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-*
Use Precision_Module
 IMPLICIT NONE
 REAL(DP), INTENT(IN)  :: R, A
 REAL(DP), INTENT(OUT) :: H (8)
 REAL(DP)              :: R2, R4

! A = PHYSICAL LENGTH OF ELEMENT     1----------2 -> R
! R = LOCAL COORDINATE OF POINT     R=0         R=1
! H = SHAPE FUNCTIONS ARRAY
! DOF ARE VALUE, SLOPE, 2ND, 3RD DERIV AT EACH END (WRT X)
! D()/DX = D()/DR DR/DX = 1/A*D()/DR

  R2 = R*R
  R4 = R2*R2
  H (1) = 1.d0 - 35.d0*R4 + 84.d0*R*R4 - 70.d0*R4*R2 + 20.d0*R4*R2*R
  H (2) = (R - 20.d0*R4 + 45.d0*R*R4 - 36.d0*R4*R2 + 10.d0*R4*R2*R)*A
  H (3) = (R2 - 10.d0*R4 + 20.d0*R*R4 - 15.d0*R4*R2 + 4.d0*R4*R2*R)*A*A/2.d0
  H (4) = (R*R2 - 4.d0*R4 + 6.d0*R*R4 - 4.d0*R4*R2 + R4*R2*R)*A*A*A/6.d0
  H (5) = 35.d0*R4 - 84.d0*R*R4 + 70.d0*R4*R2 - 20.d0*R4*R2*R
  H (6) = (10.d0*R4*R2*R - 34.d0*R4*R2 + 39.d0*R*R4 - 15.d0*R4)*A
  H (7) = (5.d0*R4 - 14.d0*R*R4 + 13.d0*R4*R2 - 4.d0*R4*R2*R)*A*A/2.d0
  H (8) = (R4*R2*R - 3.d0*R4*R2 + 3.d0*R*R4 - R4)*A*A*A/6.d0
END SUBROUTINE SHAPE_C3_L
                         
SUBROUTINE SHAPE_CUBIC (B, A, H)
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-*
! SHAPE FUNCTIONS FOR A CUBIC HERMITE, UNIT COORD
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-*
Use Precision_Module
 IMPLICIT NONE
 REAL(DP), INTENT(IN)  :: A, B
 REAL(DP), INTENT(OUT) :: H (4)

! A = LENGTH OF ELEMENT     1----------2 -> B
! B = COORDINATE OF POINT  B=0         B=1
! H = SHAPE FUNCTIONS ARRAY

  H (1) = 1.d0 - 3.d0*B*B + 2.d0*B*B*B
  H (2) = (B - 2.d0*B*B + B*B*B)*A
  H (3) = 3.d0*B*B - 2.d0*B*B*B
  H (4) = (B*B*B - B*B)*A
END SUBROUTINE SHAPE_CUBIC
                       
SUBROUTINE DERIV_SHAPE_LEGENDRE (M, X, P, DPDX)
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-*
! RECURSION RELATIONS FOR LEGENDRE POLYNOMIALS
! TO DEGREE M, AND THEIR LOCAL DERIVATIVES ON 0,1
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-*
Use Precision_Module
 IMPLICIT NONE
 INTEGER, INTENT(IN)  :: M
 REAL(DP), INTENT(IN)  :: X
 REAL(DP), INTENT(OUT) :: P (0:M), DPDX (0:M)
 INTEGER :: J
!  WARNING: THIS IS THE ONLY PLACE WHERE A ZERO SUBSCRIPT IS
!           USED IN SHAPE FUNCTIONS
  P (0) = 1.D0
  DPDX (0) = 0.D0
  IF ( M < 1) RETURN
  DPDX (1) = 1.D0
  P (1) = X
  IF ( M > 1 ) THEN
    DO J = 2, M
      P (J) = ( (J + J - 1)*X*P (J - 1) - (J - 1)*P (J - 2) )/J
      DPDX (J) = ( (J + J - 1)*X*DPDX (J - 1) - J*DPDX (J - 2)) &
               /(J - 1)
    END DO
  END IF
END SUBROUTINE DERIV_SHAPE_LEGENDRE

SUBROUTINE SHAPE_3_L (X, H)
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-*
! CALCULATE SHAPE FUNCTIONS OF A 3 NODE LINE ELEMENT
!             IN NATURAL COORDINATES
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-*
Use Precision_Module
 IMPLICIT NONE
 REAL(DP), INTENT(IN)  :: X
 REAL(DP), INTENT(OUT) :: H (3)

! H = ELEMENT SHAPE FUNCTIONS
! X = LOCAL COORDINATE OF POINT,    -1 TO +1
! LOCAL NODE COORD. ARE -1,0,+1    1-----2-----3

  H (1) = 0.5d0*(X*X - X)
  H (2) = 1.d0 - X*X
  H (3) = 0.5d0*(X*X + X)
END SUBROUTINE SHAPE_3_L
                        
SUBROUTINE SHAPE_6_T (S, T, H)
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-*
! LOCAL SHAPE FUNCTIONS FOR A SIX NODE UNIT TRIANGLE
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-*
Use Precision_Module
 IMPLICIT NONE
 REAL(DP), INTENT(IN)  :: S, T
 REAL(DP), INTENT(OUT) :: H (6)

! S,T = LOCAL COORDINATES                       3
!       OF A POINT IN THE UNIT TRIANGLE
! H = SIX SHAPE FUNCTIONS FOR A QUADRATIC  T    6 5
!     ELEMENT  WITH SIX NODES              :
! THE NODAL ORDER IS SHOWN TO THE RIGHT    .-S  1 4 2
! NODAL COORDS : 1-(0,0)   2-(1,0)   3-(0,1)
!                4-(0.5,0)  5-(0.5,0.5)  6-(0,0.5)

  H (1) = 1.d0 - 3.d0*S - 3.d0*T + 2.d0*S*S + 4.d0*S*T + 2.d0*T*T
  H (2) = 2.d0*S*S - S
  H (3) = 2.d0*T*T - T
  H (4) = 4.d0*(S - S*S - S*T)
  H (5) = 4.d0*S*T
  H (6) = 4.d0*(T - S*T - T*T)
END SUBROUTINE SHAPE_6_T
                        
SUBROUTINE SHAPE_6_T_BAR (R, S, T, H)
! * * * * * * * * * * * * * * * * * * * * * * * * * *
! BARACENTRIC SHAPES AND DERIVATIVES FOR A T6 ELEMENT
! * * * * * * * * * * * * * * * * * * * * * * * * * *
Use Precision_Module
 IMPLICIT NONE
 REAL(DP), INTENT(IN)  :: R, S, T
 REAL(DP), INTENT(OUT) :: H (6)
 REAL(DP), PARAMETER   :: ZERO = 0.D0, ONE = 1.D0, FOUR = 4.D0
                                                              
! H     = SHAPE FUNCTIONS
! R,S,T = BARACENTRIC COORDINATES, R+S+T=1
!  NODE  R  S  T       NODE   R   S   T     ORDER:
!  1     1  0  0       4     1/2 1/2  0     3
!  2     0  1  0       5      0  1/2 1/2    6 5
!  3     0  0  1       6     1/2  0  1/2    1 4 2
                                                 
!  SHAPE FUNCTIONS
  H (1) = R * (R + R - ONE)
  H (2) = S * (S + S - ONE)
  H (3) = T * (T + T - ONE)
  H (4) = FOUR * S * T
  H (5) = FOUR * R * T
  H (6) = FOUR * R * S
END SUBROUTINE SHAPE_6_T_BAR
                            
SUBROUTINE SHAPE_8_Q (S, T, H)
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-*
! SHAPE FUNCTIONS OF 8 NODE PARAMETRIC QUADRILATERAL
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-*
Use Precision_Module
 IMPLICIT NONE
 REAL(DP), INTENT(IN)  :: S, T
 REAL(DP), INTENT(OUT) :: H (8)
 REAL(DP)              :: T_plus, T_less, S_plus, S_less

! NODAL ORDER SHOWN TO RIGHT.                4 - 7 - 3
! S,T = LOCAL COORDINATES OF POINT           :   T   :
! H = SHAPE FUNCTION ARRAY                   8   *S  6
! NODE 1 AT (-1,-1)                          :       :
! NODE 3 AT (1,1)                            1 - 5 - 2

  S_plus = 1.d0 + S ; S_less = 1.d0 - S
  T_plus = 1.d0 + T ; T_less = 1.d0 - T
  H (1) = 0.25d0*S_less*T_less*(S_less + T_less - 3.d0)
  H (2) = 0.25d0*S_plus*T_less*(S_plus + T_less - 3.d0)
  H (3) = 0.25d0*S_plus*T_plus*(S_plus + T_plus - 3.d0)
  H (4) = 0.25d0*S_less*T_plus*(S_less + T_plus - 3.d0)
  H (5) = 0.5d0*T_less*(1.d0 - S*S)
  H (6) = 0.5d0*S_plus*(1.d0 - T*T)
  H (7) = 0.5d0*T_plus*(1.d0 - S*S)
  H (8) = 0.5d0*S_less*(1.d0 - T*T)
END SUBROUTINE SHAPE_8_Q
                        
SUBROUTINE SHAPE_HERMITE_RECT (N_G_DOF, HR, HS, H)
! * * * * * * * * * * * * * * * * * * * * * *
! TENSOR PRODUCT INTERPOLATION FOR HERMITES
!         ON A RECTANGULAR ELEMENT
! * * * * * * * * * * * * * * * * * * * * * *
Use Precision_Module
 IMPLICIT NONE
 INTEGER,  INTENT(IN)  :: N_G_DOF
 REAL(DP), INTENT(IN)  :: HR (2 * N_G_DOF), HS (2 * N_G_DOF)
 REAL(DP), INTENT(OUT) :: H (4 * N_G_DOF * N_G_DOF)
 INTEGER  :: I, I_PLUS, J, J_PLUS, K, KOUNT
 INTEGER  :: I_ADD (4), J_ADD (4)
 REAL(DP) :: HS_J
                                 
  DATA I_ADD / 0, 1, 1, 0 /
  DATA J_ADD / 0, 0, 1, 1 /
                           
! N_G_DOF= NO. GENERALIZED DOF PER NODE IN 1-D MODEL, >= 2
! HR, HS = 1-D HERMITES ON 2 NODES ON ORTHOGONAL SIDES
! H      = 2-D TENSOR PRODUCT FOR C^(N_G_DOF-1) CONTINUITY
!   THE FOUR  NODES ARE NUMBERED COUNTERCLOCKWISE
                                                 
  KOUNT = 0
!  LOOP OVER RECTANGLE NODES
  DO K = 1, 4
    I_PLUS = N_G_DOF * I_ADD (K)
    J_PLUS = N_G_DOF * J_ADD (K)
!    LOOP LEFT TO RIGHT IN S-DIRECTION
    DO J = 1, N_G_DOF
      HS_J = HS (J + J_PLUS)
!      LOOP LEFT TO RIGHT IN R-DIRECTION
      DO I = 1, N_G_DOF
!        FORM 2-D PRODUCT
        KOUNT     = KOUNT + 1
        H (KOUNT) = HR (I + I_PLUS) * HS_J
      END DO
    END DO
  END DO
END SUBROUTINE SHAPE_HERMITE_RECT
                   
SUBROUTINE CAUSE_VERTEX_SINGULARITY (P, H, DH, NS)
! * * * * * * * * * * * * * * * * * * * * * * * * * *
!  CONVERT STANDARD C^0 FUNCTIONS TO SINGULAR FUNCTIONS
!  WITH DERIV SINGULARITIES AT NODE NS OF O(R**(-P))
!    REF: I.J.NUM.METH.E., V10, pp 1249-1259, 1976.
! * * * * * * * * * * * * * * * * * * * * * * * * * *
!  REQUIRES SUM OF H(I) = 1, & CONST JACOBIAN
Use System_Constants
 IMPLICIT NONE
 INTEGER,  INTENT(IN)    :: NS
 REAL(DP), INTENT(IN)    :: P
 REAL(DP), INTENT(INOUT) :: H (NOD_PER_EL)
 REAL(DP), INTENT(INOUT) :: DH (N_SPACE, NOD_PER_EL)
 INTEGER  :: I, J
 REAL(DP) :: W, R
                                   
! H          = SHAPE FUNCTION ARRAY
! DH         = LOCAL DERIVATIVES OF H
! NOD_PER_EL = NO OF SHAPE FUNCTIONS
! NS         = SINGULAR NODE NUMBER,  1 <= NS <= NOD_PER_EL
! N_SPACE    = DIMENSION OF SPACE
! P          = SINGULARITY POWER, RADIUS**(-P)
                              
  IF ( P == 0.d0 ) RETURN
  W = 1.0_DP - H (NS)
  R = W**P
  DO I = 1, N_SPACE
    DO J = 1, NOD_PER_EL
      IF ( J == NS ) CYCLE ! TO NEXT NODE
      DH (I, J) = DH (I, J) / R + P * DH (I, NS) * H (J) / R / W
    END DO
    DH (I, NS) = (1.0_DP - P) * DH (I, NS) / R
  END DO
  DO J = 1, NOD_PER_EL
    IF (J == NS ) CYCLE ! TO NEXT NODE
    H (J) = H (J) / R
  END DO
  H (NS) = 1.0_DP - W / R
END SUBROUTINE CAUSE_VERTEX_SINGULARITY
                     
SUBROUTINE CAUSE_NODE_SINGULARITY (P, H, DH, NS) 
! * * * * * * * * * * * * * * * * * * * * * * * * * *               
!  CONVERT STANDARD C^0 FUNCTIONS TO SINGULAR FUNCTIONS             
!  WITH DERIV SINGULARITIES AT NODE NS OF O(R**(-P))                
!    REF: I.J.NUM.METH.E., V10, pp 1249-1259, 1976.
! * * * * * * * * * * * * * * * * * * * * * * * * * *               
Use System_Constants ! for N_SPACE, NOD_PER_EL
IMPLICIT NONE
  INTEGER,  INTENT(IN)    :: NS
  REAL(DP), INTENT(IN)    :: P
  REAL(DP), INTENT(INOUT) :: H (NOD_PER_EL)
  REAL(DP), INTENT(INOUT) :: DH (N_SPACE, NOD_PER_EL) 
  INTEGER  :: I, J
  REAL(DP) :: R, W

! H          = SHAPE FUNCTION ARRAY
! DH         = LOCAL DERIVATIVES OF H                 
! NOD_PER_EL = NO OF SHAPE FUNCTIONS
! N_SPACE    = DIMENSION OF SPACE                
! NS         = SINGULAR NODE NUMBER,  1 <= NS <= NOD_PER_EL      
! P          = SINGULARITY POWER, RADIUS**(-P)

!   REQUIRES SUM OF H(I) = 1, & CONST JACOBIAN                        
   IF (P == 0.0_DP) RETURN 
     W = 1.0_DP - H (NS) 
     R = W**P 
     DO I = 1, N_SPACE 
       DO J = 1, NOD_PER_EL 
         IF (J == NS) CYCLE ! J
           DH (I, J) = DH (I, J) / R &
                     + P * DH (I, NS) * H (J) / R / W 
       END DO 
       DH (I, NS) = (1.0_DP - P) * DH (I, NS) / R 
     END DO 
     DO J = 1, NOD_PER_EL 
       IF (J == NS) CYCLE ! J
         H (J) = H (J) / R 
     END DO 
     H (NS) = 1.0_DP - W / R 
END SUBROUTINE CAUSE_NODE_SINGULARITY                         

FUNCTION GET_SCALAR_DEGREE ()  RESULT (DEGREE_L)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
! RETURN THE "COMPLETE" DEGREE IN AN ELEMENT INTERPOLATION
!            (FOR USE BY ERROR ESTIMATORS)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants
Use Elem_Type_Data
 IMPLICIT NONE
 INTEGER            :: DEGREE_L   
 INTEGER, PARAMETER :: ONE = 1    
                                                  
! DEGREE_L = COMPLETE (NOT EDGE) DEGREE OF AN ELEMENT
! LT_N     = NUMBER OF NODES PER ELEMENT
! LT_PARM  = NUMBER OF SPATIAL DIMENSIONS
! N_G_DOF  = NUMBER OF DEGREES OF FREEDOM PER NODE

   DEGREE_L = 1

!  BRANCH ON SPACE, THEN NUMBER OF NODES 
  SELECT CASE ( LT_PARM )  ! LOCAL SPACE DIMENSION
    CASE (0) ! point
      PRINT *,'WARNING, NO POINT ELEMENT DEGREE'
      N_WARN = N_WARN + 1 ; DEGREE_L = 0
    CASE (1) ! 1-D space

!-->     1-D ELEMENTS
      SELECT CASE ( LT_N ) ! NUMBER OF ELEMENT NODES
        CASE (2) ; DEGREE_L = 1
        CASE (3) ; DEGREE_L = 2
        CASE (4) ; DEGREE_L = 3
        CASE DEFAULT 
          CALL DEGREE_USER (DEGREE_L, LT_N, LT_PARM, ONE)     
      END SELECT ! number of element nodes
    CASE (2) ! 2-D space
      SELECT CASE ( LT_N )

!-->       TRIANGULAR 2-D ELEMENTS
        CASE ( 3) ; DEGREE_L = 1
!       CASE ( 4) ; DEGREE_L = 1
        CASE ( 6) ; DEGREE_L = 2
!       CASE ( 7) ; DEGREE_L = 2
        CASE (10) ; DEGREE_L = 3
        CASE (15) ; DEGREE_L = 4

!-->       QUADRILATERAL 2-D ELEMENTS
        CASE ( 4) ; DEGREE_L = 1
        CASE ( 8) ; DEGREE_L = 2
        CASE ( 9) ; DEGREE_L = 2
        CASE (12) ; DEGREE_L = 3 
        CASE (16) ; DEGREE_L = 3 ! IF LAGRANGIAN, ELSE = 2
!       CASE (17) ; DEGREE_L = 4
!       CASE (25) ; DEGREE_L = 4
        CASE DEFAULT 
          CALL DEGREE_USER (DEGREE_L, LT_N, LT_PARM, ONE)     
      END SELECT ! number of element nodes
    CASE (3) ! 3-D space
      SELECT CASE ( LT_N )

!-->       HEXAHEDRA 3-D ELEMENTS
        CASE ( 8) ; DEGREE_L = 1 
        CASE (20) ; DEGREE_L = 2
!       CASE (27) ; DEGREE_L = 3
!       CASE (32) ; DEGREE_L = 3

!-->       TETRAHEDRA 3-D ELEMENTS (PYRAMIDS)
!       CASE ( 4) ; DEGREE_L = 1
!       CASE (10) ; DEGREE_L = 2
!       CASE (20) ; DEGREE_L = 3

!-->       WEDGE 3-D ELEMENTS
!       CASE ( 6) ; DEGREE_L = 1
!       CASE (15) ; DEGREE_L = 2

!-->       USER SUPPLIED ELEMENT
        CASE DEFAULT 
          CALL DEGREE_USER (DEGREE_L, LT_N, LT_PARM, ONE)     
      END SELECT ! number of element nodes
    CASE DEFAULT
      PRINT *, LT_N, LT_PARM, N_G_DOF
      STOP 'INVALID SPACE SIZE IN GET_SCALAR_DEGREE'
END SELECT ! number of spatial dimensions
END FUNCTION GET_SCALAR_DEGREE

SUBROUTINE DEGREE_USER (DEGREE_L, LT_N, LT_PARM, N_G_DOF)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
! RETURN THE "COMPLETE" DEGREE IN AN ELEMENT INTERPOLATION
!            (FOR USE BY ERROR ESTIMATORS)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
 IMPLICIT NONE
 INTEGER, INTENT(IN)  :: LT_N, LT_PARM, N_G_DOF
 INTEGER, INTENT(OUT) :: DEGREE_L   
                                                  
! DEGREE_L = COMPLETE (NOT EDGE) DEGREE OF AN ELEMENT
! LT_N     = NUMBER OF NODES PER ELEMENT
! LT_PARM  = NUMBER OF SPATIAL DIMENSIONS
! N_G_DOF  = NUMBER OF DEGREES OF FREEDOM PER NODE
  DEGREE_L = 0
  PRINT *,'DEGREE_USER: LT_N, LT_PARM, N_G_DOF ', LT_N, LT_PARM, N_G_DOF
  STOP 'NO USER DEFINITION IN DEGREE_USER'
END SUBROUTINE DEGREE_USER 

FUNCTION GET_C_N_DEGREE ()  RESULT (DEGREE_L)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
! RETURN THE "COMPLETE" DEGREE IN A C_N ELEMENT INTERPOLATION
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants
Use Elem_Type_Data
 IMPLICIT NONE
 INTEGER :: DEGREE_L     	
                                                  
! DEGREE_L = COMPLETE (NOT EDGE) DEGREE OF AN ELEMENT	
! LT_N     = NUMBER OF NODES PER ELEMENT
! LT_PARM  = NO OF SPATIAL DIMENSIONS
! N_G_DOF  = NUMBER OF DEGREES OF FREEDOM PER NODE

!  BRANCH ON SPACE, THEN NUMBER OF NODES 
  SELECT CASE ( LT_PARM )
    CASE (0) ! point
      PRINT *,'WARNING, NO POINT ELEMENT C_N DEGREE'
      N_WARN = N_WARN + 1 ; DEGREE_L = 0
    CASE (1) ! 1-D space

!-->     1-D ELEMENTS
      SELECT CASE ( LT_N )
        CASE (2) 
          SELECT CASE ( N_G_DOF )
            CASE (2) ; DEGREE_L = 3   ! FOR SHAPE_C1_L 
            CASE (3) ; DEGREE_L = 5   ! FOR SHAPE_C2_L 
            CASE (4) ; DEGREE_L = 7   ! FOR SHAPE_C3_L 
            CASE DEFAULT 
              CALL DEGREE_C_N_USER (DEGREE_L, LT_N, &
                                  LT_PARM, N_G_DOF)
          END SELECT
        CASE DEFAULT 
          CALL DEGREE_C_N_USER (DEGREE_L, LT_N, &
                              LT_PARM, N_G_DOF)
      END SELECT ! number of element nodes
    CASE (2) ! 2-D space
      SELECT CASE ( LT_N )

!-->       RECTANGULAR 2-D ELEMENTS
!       CASE (16) 
!         DEGREE_L = 3   ! FOR SHAPE_16_R 
        CASE DEFAULT 
          CALL DEGREE_C_N_USER (DEGREE_L, LT_N, LT_PARM, N_G_DOF)
      END SELECT ! number of element nodes

    CASE DEFAULT  ! 3D, OR ERROR
      PRINT *, LT_N, LT_PARM, N_G_DOF
      STOP 'INVALID DEGREE IN GET_C_N_DEGREE'
END SELECT ! number of spatial dimensions
END FUNCTION GET_C_N_DEGREE

SUBROUTINE DEGREE_C_N_USER (DEGREE_L, LT_N, LT_PARM, N_G_DOF)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
! RETURN THE "COMPLETE" DEGREE IN A C_1 ELEMENT INTERPOLATION
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
 IMPLICIT NONE
 INTEGER, INTENT(IN)  :: LT_N, LT_PARM, N_G_DOF
 INTEGER, INTENT(OUT) :: DEGREE_L     	
                                                  
! DEGREE_L = COMPLETE (NOT EDGE) DEGREE OF AN ELEMENT	
! LT_N     = NUMBER OF NODES PER ELEMENT
! LT_PARM  = NO OF SPATIAL DIMENSIONS
! N_G_DOF  = NUMBER OF DEGREES OF FREEDOM PER NODE
  DEGREE_L = 0
  PRINT *,'DEGREE_C_N_USER:', LT_N, LT_PARM, N_G_DOF
  STOP 'NO USER DEFINITION IN DEGREE_C_N_USER'
END SUBROUTINE DEGREE_C_N_USER 

SUBROUTINE SCALAR_DEGREES (DEGREE)   
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
! RETURN DEGREE OF THE COMPLETE POLYNOMIAL COMTAINED IN ANY
!       C_0 ELEMENT INTERPOLATION FUNCTION
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants
Use Elem_Type_Data
 IMPLICIT NONE
 INTEGER, INTENT(OUT) :: DEGREE
                                                  
! DEGREE  = COMPLETE POLYNOMIAL ORDER
! LT_N    = NUMBER OF NODES PER ELEMENT
! LT_PARM = NO OF SPATIAL DIMENSIONS
! N_G_DOF = NUMBER OF DEGREES OF FREEDOM PER NODE

!  BRANCH ON SPACE, THEN NUMBER OF NODES 
  SELECT CASE ( LT_PARM )

    CASE (0) ! point
      PRINT *,'WARNING, NO POINT ELEMENT C_0 DEGREE'
      N_WARN = N_WARN + 1 ; DEGREE = 0
    CASE (1) ! 1-D space
!-->     1-D ELEMENTS
      SELECT CASE ( LT_N )
        CASE (2) 
          SELECT CASE ( N_G_DOF )
            CASE (1) ; DEGREE = 1   ! SHAPE_2_L
            CASE DEFAULT 
              CALL DEGREE_USER (DEGREE, LT_N, LT_PARM, N_G_DOF)
          END SELECT
        CASE (3) ; DEGREE = 2   ! SHAPE_3_L
        CASE (4) ; DEGREE = 3   ! SHAPE_4_L
        CASE DEFAULT 
          CALL DEGREE_USER (DEGREE, LT_N, LT_PARM, N_G_DOF)
      END SELECT ! number of element nodes

    CASE (2) ! 2-D space
      SELECT CASE ( LT_N )
!-->       TRIANGULAR 2-D ELEMENTS
        CASE ( 3) ; DEGREE = 1   ! SHAPE_3_T
!       CASE ( 4) ; DEGREE = 1   ! SHAPE_4_T
        CASE ( 6) ; DEGREE = 2   ! SHAPE_6_T
!       CASE ( 7) ; DEGREE = 2   ! SHAPE_7_T
        CASE (10) ; DEGREE = 3   ! SHAPE_10_T
        CASE (15) ; DEGREE = 4   ! SHAPE_15_T

!-->       QUADRILATERAL 2-D ELEMENTS
        CASE ( 4) ; DEGREE = 1   ! SHAPE_4_Q
        CASE ( 8) ; DEGREE = 2   ! SHAPE_8_Q
        CASE ( 9) ; DEGREE = 2   ! SHAPE_9_Q
        CASE (12) ; DEGREE = 3   ! SHAPE_12_Q
        CASE (16) ; DEGREE = 3   ! SHAPE_16_Q, LAGRANGIAN
        CASE (17) ; DEGREE = 4   ! SHAPE_17_Q
!       CASE (25) ; DEGREE = 4   ! SHAPE_25_Q
        CASE DEFAULT 
          CALL DEGREE_USER (DEGREE, LT_N, LT_PARM, N_G_DOF)
      END SELECT ! number of element nodes

    CASE (3) ! 3-D space
      SELECT CASE ( LT_N )
!-->       HEXAHEDRA 3-D ELEMENTS
        CASE ( 8) ; DEGREE = 1   ! SHAPE_8_H 
        CASE (20) ; DEGREE = 2   ! SHAPE_20_H
!       CASE (27) ; DEGREE = 3   ! SHAPE_27_H
!       CASE (32) ; DEGREE = 3   ! SHAPE_32_H

!-->       TETRAHEDRA 3-D ELEMENTS (PYRAMIDS)
        CASE ( 4) ; DEGREE = 1   ! SHAPE_4_P
!       CASE (10) ; DEGREE = 2   ! SHAPE_10_P
!       CASE (20) ; DEGREE = 3   ! SHAPE_20_P

!-->       WEDGE 3-D ELEMENTS
!       CASE ( 6) ; DEGREE = 1   ! SHAPE_6_W
!       CASE (15) ; DEGREE = 2   ! SHAPE_15_W

!-->       USER SUPPLIED ELEMENT
        CASE DEFAULT 
          CALL DEGREE_USER (DEGREE, LT_N, LT_PARM, N_G_DOF)
      END SELECT ! number of element nodes

    CASE DEFAULT 
      PRINT *, LT_N, LT_PARM, N_G_DOF
      STOP 'INVALID SPACE SIZE IN SCALAR_DEGREES'
 END SELECT ! number of spatial dimensions
END SUBROUTINE SCALAR_DEGREES
                         
SUBROUTINE SHAPE_C1_L3_NC (P, A, H)
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-*
!                    untested 8/2004
! SHAPE FUNCTIONS FOR QUINTIC HERMITE IN NATURAL COORDINATES
! ( A 3 NODED C1 ELEMENT, REF: L. R. Ram-Mohn 2002 OUP )
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-*
Use Precision_Module
 IMPLICIT NONE
 REAL(DP), INTENT(IN)  :: P, A
 REAL(DP), INTENT(OUT) :: H (6)
 REAL(DP) :: P_2, P_3, P_4, P_5

! A = PHYSICAL LENGTH OF ELEMENT     1------2------3 ---> P
! P = LOCAL COORDINATE OF POINT    P=-1    P=0     P=1
! H = SHAPE FUNCTIONS ARRAY
! DOF ARE FUNCTION AND SLOPE, WRT X, AT EACH NODE (N_G_DOF=2)
! D()/DX = D()/DP DP/DX = 2/A * D()/DP

  P_2 = P * P   ; P_3 = P * P_2
  P_4 = P * P_3 ; P_5 = P * P_4
  H(1) = (4*P_2 - 5*P_3 - 2*P_4 + 3*P_5) * 0.25d0
  H(2) = (P_2 - P_3 - P_4 + P_5) * 0.125d0 * A
  H(3) =  1 - 2*P_2 + P_4
  H(4) = (P - 2*P_3 + P_5) * A * 0.5d0
  H(5) = (4*P_2 + 5*P_3 - 2*P_4 - 3*P_5) * 0.25d0 
  H(6) = (-P_2 - P_3 + P_4 + P_5) * 0.125d0 * A
END SUBROUTINE SHAPE_C1_L3_NC

SUBROUTINE DERIV_C1_L3_NC (P, A, DH)
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-*
!                    untested 8/2004
!      FIRST PHYSICAL DERIVATIVES OF QUINTIC HERMITE 
!               IN NATURAL COORDINATES
! ( A 3 NODED C1 ELEMENT, REF: L. R. Ram-Mohn 2002 OUP )
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-*
Use Precision_Module
 IMPLICIT NONE
 REAL(DP), INTENT(IN)  :: P, A
 REAL(DP), INTENT(OUT) :: DH (6)
 REAL(DP) :: P_2, P_3, P_4

! A = PHYSICAL LENGTH OF ELEMENT        1------2------3 ---> P 
! R = LOCAL COORDINATE OF POINT        P=-1    P=0     P=1
! DH = FIRST PHYSICAL DERIVATIVES OF H
! D()/DX = D()/DP DP/DX = 2/A * D()/DP
! DOF ARE FUNCTION AND SLOPE, D/DX, AT EACH NODE

  P_2 = P * P   ; P_3 = P * P_2 ; P_4 = P * P_3
  DH (1) = (8*P - 15*P_2 - 8*P_3 + 15*P_4) * 0.5d0 / A
  DH (2) = (2*P - 3*P_2 - 4*P_3 + 5*P_4) * 0.25d0
  DH (3) = (-4*P + 4*P_3) * 2 / A
  DH (4) = (1 - 6*P_2 + 5*P_4) 
  DH (5) = (8*P + 15*P_2 - 8*P_3 - 15*P_4) * 0.5d0 / A
  DH (6) = (-2*P - 3*P_2 + 4*P_3 + 5*P_4) * 0.25d0
END SUBROUTINE DERIV_C1_L3_NC

SUBROUTINE DERIV2_C1_L3_NC (P, A, D2H)
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-*
!                    untested 8/2004
!      SECOND PHYSICAL DERIVATIVES OF QUINTIC HERMITE 
!               IN NATURAL COORDINATES
! ( A 3 NODED C1 ELEMENT, REF: L. R. Ram-Mohn 2002 OUP )
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-*
Use Precision_Module
 IMPLICIT NONE
 REAL(DP), INTENT(IN)  :: P, A
 REAL(DP), INTENT(OUT) :: D2H (6)
 REAL(DP) :: P_2, P_3

! A = PHYSICAL LENGTH OF ELEMENT        1------2------3 ---> P 
! R = LOCAL COORDINATE OF POINT        P=-1    P=0     P=1
! D2H = SECOND PHYSICAL DERIVATIVES OF H
! D()/DX = D()/DP DP/DX = 2/A * D()/DP
! DOF ARE FUNCTION AND SLOPE, D/DX, AT EACH NODE

  P_2 = P * P   ; P_3 = P * P_2
  D2H (1) = (8 - 30*P - 24*P_2 + 60*P_3) / A **2
  D2H (2) = (2 - 6*P - 12*P_2 + 20*P_3) * 0.5d0 / A
  D2H (3) = (-4 + 12*P_2) * 4 / A **2
  D2H (4) = (-12*P + 20*P_3) * 2 / A 
  D2H (5) = (8 + 30*P - 24*P_2 - 60*P_3) / A **2
  D2H (6) = (-2 - 6*P + 12*P_2 + 20*P_3) * 0.5d0 / A
END SUBROUTINE DERIV2_C1_L3_NC
