! copyright 2005,  J. E. Akin, all rights reserved.
SUBROUTINE BARACENTRIC_TO_UNIT_DERIV (N_GEOM, N_PARM, DLB, DLU)
! * * * * * * * * * * * * * * * * * * * * * * * * * * *
!    CONVERT SIMPLEX BARACENTRIC LOCAL DERIVATIVES TO
!        UNIT COORDINATE LOCAL DERIVATIVES
! * * * * * * * * * * * * * * * * * * * * * * * * * * *
Use Precision_Module ! for DP
 IMPLICIT NONE
 INTEGER,  INTENT(IN) :: N_GEOM, N_PARM
 REAL(DP), INTENT(IN) :: DLB (N_PARM+1, N_GEOM)
 REAL(DP), INTENT(OUT):: DLU (N_PARM, N_GEOM)
 INTEGER  :: J, K
 REAL(DP) :: DLB_LAST
                                    
! DLB    = LOCAL DERIVATIVES IN BARACENTRIC COORDINATES
! DLU    = LOCAL DERIVATIVES IN UNIT COORDINATES
! N_GEOM = NUMBER OF NODES
! N_PARM = DIMENSION OF UNIT SIMPLEX SPACE
                                           
!  LOOP OVER THE INTERPOLATION FUNCTIONS
  DO K = 1, N_GEOM
    DLB_LAST = DLB (N_PARM+1, K)
!    SUBTRACT OFF THE LAST BARACENTRIC VALUE
    DO J = 1, N_PARM
      DLU (J, K) = DLB (J, K) - DLB_LAST
    END DO
  END DO
END SUBROUTINE BARACENTRIC_TO_UNIT_DERIV
                                  
SUBROUTINE FUND_MAG_ONE (N_GEOM, DELTA, COORD, RJ, A, FFM)
! * * * * * * * * * * * * * * * * * * * * * * * * * *
! FIND 1ST FUNDAMENTAL MAGNITUDES ON A 3-D SURFACE
! * * * * * * * * * * * * * * * * * * * * * * * * * *
Use Precision_Module
 IMPLICIT NONE
 INTEGER,  INTENT(IN)  :: N_GEOM
 REAL(DP), INTENT(IN)  :: DELTA (2, N_GEOM)
 REAL(DP), INTENT(IN)  :: COORD (N_GEOM, 3)
 REAL(DP), INTENT(OUT) :: A (2, 2), RJ (2, 3), FFM
                                                  
! N_GEOM = NUMBER OF NODES DEFINING THE SURFACE
! DELTA  = LOCAL DERIVATIVES OF THE SHAPE FUNCTIONS
! COORD  = THREE SPATIAL COORDINATES OF THE NODES
! RJ     = REDUCED JACOBIAN MATRIX
! A      = 1ST FUND. MAG. (METRIC TENSOR)
! FFM    = DET(A)

! INDEX NOTATION:  I=1,2,3,  R,S=1,2,  _=SUBSCRIPTS
!             A_RS = X_I,R * X_I,S  RJ_RI = X_I,R

!  FORM REDUCED JACOBIAN
  RJ = MATMUL (DELTA, COORD)

!  FORM METRIC TENSOR
  A (1, 1) = RJ (1, 1) * RJ (1, 1) + RJ (1, 2) * RJ (1, 2) &
           + RJ (1, 3) * RJ (1, 3)
  A (1, 2) = RJ (1, 1) * RJ (2, 1) + RJ (1, 2) * RJ (2, 2) &
           + RJ (1, 3) * RJ (2, 3)
  A (2, 2) = RJ (2, 1) * RJ (2, 1) + RJ (2, 2) * RJ (2, 2) &
           + RJ (2, 3) * RJ (2, 3)
  FFM = A (1, 1) * A (2, 2) - A (1, 2) * A (1, 2)

! SURFACE AREA DA = SQRT(FFM)*DR*DS
! ANGLE R-S    COS (ANG) = A12/SQRT( A11*A22 )
! LENGTHS      DL1 = SQRT(A11)*DR ;  DL2 = SQRT(A22)*DS
  IF ( FFM < 0.0 ) PRINT *,'FATAL DATA ERROR IN FUND_MAG_ONE'
END SUBROUTINE FUND_MAG_ONE
                           
FUNCTION INSIDE_2D_POLYGON (N_VERT, XY, COORD) RESULT (IN_CODE)
!* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
! FIND IF PT (X,Y) IS INSIDE ANY 2-D POLYGON WITH N_VERT VERTICIES
!* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
!          IEEE TRAN EL COMP, EC-16, NO-6, DEC, 67
 Use Precision_Module
  IMPLICIT NONE
  INTEGER,  INTENT(IN) :: N_VERT             ! number of vertices
  REAL(DP), INTENT(IN) :: XY (2)             ! coordinates of point
  REAL(DP), INTENT(IN) :: COORD (N_VERT, 2)  ! coords of vertices
  INTEGER              :: IN_CODE            ! return status
  REAL(DP), PARAMETER  :: TWO_PI = 6.28318530717958623_DP
  REAL(DP), PARAMETER  :: ZERO   = 2*TINY (1.0_DP)
  REAL(DP) :: ALFA (N_VERT + 1), BETA (N_VERT + 1), DIFF ! work 
  INTEGER  :: I

! FOR RIGHT HAND SYSTEM ORDER VERTICIES IN CCW MANNER
! IN_CODE = 0 IF OUTSIDE, IN_CODE = 1 IF INSIDE, = ?? if on edge
   IN_CODE = 0
   DO I = 1, N_VERT
     ALFA (I) = ATAN2 ((COORD (I, 2) - XY (2)), &
                       (COORD (I, 1) - XY (1)))/TWO_PI
   END DO 
   ALFA (N_VERT + 1) = ALFA (1) ; BETA (1)     = ALFA (1)
   DO I = 2, N_VERT + 1
     DIFF = ALFA (I) - ALFA (I-1)
     IF ( ABS (DIFF - 0.5_DP) <= ZERO )  THEN    
       IN_CODE = 1 ; RETURN
     ELSE
       BETA (I) = BETA (I-1) + DIFF
       IF ( DIFF >   0.5_DP )  BETA (I) = BETA (I) - 1.0_DP
       IF ( DIFF < (-0.5_DP) ) BETA (I) = BETA (I) + 1.0_DP
     END IF
   END DO  
   IF ( ABS (BETA (N_VERT + 1) - (BETA (1) + 1.0_DP)) <= ZERO ) &
     IN_CODE = 1
END FUNCTION INSIDE_2D_POLYGON 

SUBROUTINE INSIDE_3T (XY, COORD, INSIDE)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
!  IS POINT XY INSIDE A TRIANGLE WITH GIVEN COORDINATES
!  (NODES ARE ASSUMED TO BE NUMBERED COUNTERCLOCKWISE)
!            (SEE INSIDE_2D_POLYGON AS ALTERNATE)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
Use Precision_Module
 IMPLICIT NONE
 REAL(DP), INTENT(IN)  :: XY (2), COORD (3, 2)
 INTEGER,  INTENT(OUT) :: INSIDE
 REAL(DP) :: TWO_A
                                  
! COORD  = PHYSICAL COORDINATES OF NODES
! INSIDE = 0 IF XY IS NOT IN ELEMENT, ELSE = 1
! XY     = PHSYICAL COORDINATES OF POINT

  INSIDE = 0 !  INITALIZE AS OUTSIDE (ON EDGE IS INSIDE)

!  TESTING TRIANGULAR AREA ON 3 SIDES (+ IF IN,0 ON EDGE)
  TWO_A = COORD (1, 1) * (COORD (2, 2) - XY (2) ) &
        + COORD (2, 1) * (XY (2) - COORD (1, 2) ) &
        + XY (1) * (COORD (1, 2) - COORD (2, 2) )
  IF ( TWO_A < 0.0 ) RETURN
    TWO_A = COORD (2, 1) * (COORD (3, 2) - XY (2) ) &
          + COORD (3, 1) * (XY (2) - COORD (2, 2) ) &
          + XY (1) * (COORD (2, 2) - COORD (3, 2) )
    IF ( TWO_A < 0.0 ) RETURN
      TWO_A = COORD (3, 1) * (COORD (1, 2) - XY (2) ) &
            + COORD (1, 1) * (XY (2) - COORD (3, 2) ) &
            + XY (1) * (COORD (3, 2) - COORD (1, 2) )
      IF ( TWO_A < 0.0 ) RETURN
!          POINT IS INSIDE
          INSIDE = 1
END SUBROUTINE INSIDE_3T
                        
SUBROUTINE IS_INSIDE (XYZ, INSIDE)
! * * * * * * * * * * * * * * * * * * * * * * * * *
! is xyz inside a line, or flat T3 or Q4 element ?
! * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants
Use Elem_Type_Data  ! for GEOMETRY (LT_GEOM, N_SPACE)
  REAL(DP), INTENT(IN)  :: XYZ (N_SPACE) ! point COORD
  INTEGER,  INTENT(OUT) :: INSIDE  ! 0=outside, 1=inside or on edge
  REAL(DP) :: Q_134 (3, n_space)   ! 2nd triangle of quadrilateral
            
! GEOMETRY = COORDINATES OF ELEMENT GEOMETRY NODES
! INSIDE   = FLAG FOR POINT LOCATION
! LT_GEOM  = NUMBER OF GEOMETRICAL NODES
! N_SPACE  = DIMENSION OF PHYSICAL SPACE
! XYZ      = COORDINATES OF TEST POINT

    inside = 0                                  ! outside element ?
    select case (N_SPACE)
    case (1)
      if ( xyz(1) > minval( GEOMETRY (:,1) ) .and. &
           xyz(1) < maxval( GEOMETRY (:,1) ) ) inside = 1
    case (2)
      ! see INSIDE_2D_POLYGON as alternative
      call INSIDE_3T (xyz, GEOMETRY(1:3,:), inside)   ! check triangle
      if ( inside == 1 ) return                    ! was there
      if ( LT_N == 4 .or. LT_N == 8 .or. &    ! check quadrilateral
           LT_N == 9 ) then
        ! also see WHERE_IN_4Q for alternative
        q_134 (1,:) = GEOMETRY (1,:)             ! build other triangle
        q_134 (2,:) = GEOMETRY (3,:)
        q_134 (3,:) = GEOMETRY (4,:)
        call INSIDE_3T (xyz, q_134, inside)   ! check other triangle
      end if
    case default
      stop 'Extend is_inside to 3-D'
    end select
END SUBROUTINE IS_INSIDE
                          
SUBROUTINE LCONTR_ (R_IN_OUT, S_IN_OUT, V, COORD, XYZ, N_GEOM, &
                    N_SPACE, KOUNT, L_TYPE)
!  * * * * * * * * * * * * * * * * * * * * * * * * * *
!  CALCULATE COORDINATES OF POINTS ON A CONTOUR
!    ON AN ISOPARAMETRIC SURFACE
!  * * * * * * * * * * * * * * * * * * * * * * * * * *
Use Precision_Module
 IMPLICIT NONE
 INTEGER,  PARAMETER     :: MAX = 300
 REAL(DP), PARAMETER     :: TRY_DL = 0.005
 INTEGER,  INTENT(IN)    :: N_GEOM, N_SPACE, L_TYPE
 REAL(DP), INTENT(IN)    :: V (N_GEOM)
 REAL(DP), INTENT(IN)    :: COORD (N_GEOM, N_SPACE)
 INTEGER,  INTENT(OUT)   :: KOUNT
 REAL(DP), INTENT(OUT)   :: XYZ (MAX, N_SPACE)
 REAL(DP), INTENT(INOUT) :: R_IN_OUT, S_IN_OUT
                                              
 INTEGER  :: I
 REAL(DP) :: R (MAX), S (MAX), DER (2), PT(N_SPACE), P_XYZ (3)
 REAL(DP) :: H (N_GEOM), DH (2, N_GEOM)
 REAL(DP) :: DL, DV_DR, DV_DS, GRAD, R_NEW, S_NEW
                                                          
 EQUIVALENCE (DER (1), DV_DR), (DER (2), DV_DS)
                                                
!  MAX      = MAX. NO. CONTOUR POINTS
!  DL       = LOCAL COORD. LENGTH OF EACH SEGMENT
!  KOUNT    = NO. OF PTS. ON CONTOUR CURVE
!  R,S      = LOCAL COORDINATES
!  XYZ      = GLOBAL COORD. ARRAY FOR COUNTOUR PTS.
!  PT       = LOCAL  COORDINATES OF A SINGLE POINT
!  P_XYZ    = GLOBAL COORDINATES OF A SINGLE POINT
!  COORD    = GLOBAL COORD. OF NODES OF ELEMENT
!  N_GEOM   = NUMBER OF NODES PER ELEMENT
!  N_SPACE  = DIMENSION OF GLOBAL SPACE, 2 OR 3
!  R_IN_OUT = R COORD OF 1ST (IN) OR LAST (OUT) PT
!  S_IN_OUT = S COORD OF 1ST (IN) OR LAST (OUT) PT
!  H        = ELEMENT INTERPOLATION FUNCTIONS
!  DH       = LOCAL COORD. DERIVATIVES OF H
!  V        = NODAL VALUES OF QUANTITY TO BE CONTOURED
!  L_TYPE   = ELEM. TYPE, 0=QUADRILATERAL, 1=TRIANGLE
                                                    
  KOUNT = 1
  DL    = TRY_DL
  IF ( L_TYPE == 0 ) DL = DL * 2.0

!-->    LOCAL COORDINATE CALCULATIONS
  R (1) = R_IN_OUT ; S (1) = S_IN_OUT

!  MARCH ALONG CONTOUR
  5 CONTINUE  !9 use do while loop

!  FORM SHAPE FUNCTION LOCAL DERIVATIVES
  PT (1) = R (KOUNT) ; PT (2) = S (KOUNT)
  CALL SCALAR_DERIVS (PT, DH)

!  FORM LOCAL DERIVATIVE OF VARIABLE, DER = DH*V
  DER = MATMUL ( DH, V ) 
  GRAD = SQRT (DV_DR * DV_DR + DV_DS * DV_DS)

!  FIND LOCAL COORD. OF SEGMENT END
  R_NEW = R (KOUNT) - DL * DV_DS / GRAD
  S_NEW = S (KOUNT) + DL * DV_DR / GRAD

!  IS NEXT POINT IN THE ELEMENT
  IF ( L_TYPE == 0 ) THEN
!     QUADRILATERAL (-1 TO +1): IS POINT OUTSIDE ?
10   IF ( ABS (R_NEW) > 1.0 .OR. ABS (S_NEW) > 1.0 ) GO TO 20
  ELSE
!     TRIANGLE (UNIT COORDINATES): OUTSIDE ?
    IF ( R_NEW < 0.0 .OR. S_NEW < 0.0 .OR. &
        (R_NEW + S_NEW) > 1.0) GO TO  20
  END IF

!  ADD POINT TO CONTOUR LIST
  KOUNT     = KOUNT + 1
  R (KOUNT) = R_NEW
  S (KOUNT) = S_NEW
  IF ( KOUNT < MAX ) GO TO 5

!-->    LINE COMPLETED IN LOCAL COORDINATES
20 CONTINUE  !9 use do while loop
  IF ( KOUNT == 1 ) THEN
    DL = - DL
    GO TO 5
  END IF

!  RETURN LOCAL COORD. OF LAST POINT
  R_IN_OUT = R (KOUNT)
  S_IN_OUT = S (KOUNT)

!-->    CONVERT TO GLOBAL COORDINATES
  DO I = 1, KOUNT

!     EVALUATE SHAPE FUNCTIONS
    CALL SCALAR_SHAPES (PT, H)

!     EVALUATE GLOBAL COORDINATES, XYZ = H*COORD
    P_XYZ = MATMUL ( H, COORD ) 
    XYZ (I, 1) = P_XYZ (1)
    XYZ (I, 2) = P_XYZ (2)
    IF ( N_SPACE == 3 ) XYZ (I, 3) = P_XYZ (3)
  END DO
!-->    CONTOUR FINISHED, RETURN FOR PLOTTING
END SUBROUTINE LCONTR_
                          
SUBROUTINE SURFACE_LAME (GRAD, ROOT_E, F, ROOT_G, H)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
!     GET THE SURFACE LAME PARAMETERS, E F G AND H,
!     AT A POINT ON A PARAMETRIC CURVE OR SURFACE
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants
 IMPLICIT NONE
 REAL(DP), INTENT(IN)  :: GRAD (N_PARM, N_GEOM)
 REAL(DP), INTENT(OUT) :: ROOT_E, F, ROOT_G, H
 REAL(DP) :: E, G
                                 
! NOD_PER_EL = NUMBER OF NODES PER ELEMENT
! N_PARM     = NUMBER OF PARAMETRIC DIMENSIONS <= N_GEOM
! N_GEOM     = DIMENSION OF PHYSICAL SPACE
! DELTA      = LOCAL DERIVATIVES OF NOD_PER_EL INTERPOLATION FUNCTIONS 
! COORD      = SPATIAL COORDINATES OF ELEMENT'S NODES
! GRAD       = TANGENT MATRIX FROM TANGENT_VECTOR; 
!              = JACOBIAN IF N_PARM = N_GEOM
! ROOT_E     = SQUARE ROOT OF E, E = X,R^2   + Y,R^2   + Z,R^2
! ROOT_F     = SQUARE ROOT OF F, F = X,S*X,R + Y,S*Y,R + Z,S*Z,R
! ROOT_G     = SQUARE ROOT OF G, G = X,S^2   + Y,S^2   + Z,S^2
                                                          
  IF ( N_PARM > 2) STOP 'INVALID el_parm ARGUMENT, SURFACE_LAME'
  G = 0. ; F = 0. ;  E =     GRAD (1, 1)**2
  IF ( N_GEOM >  1) E = E + GRAD (1, 2)**2
  IF ( N_GEOM == 3) E = E + GRAD (1, 3)**2
  IF ( N_PARM > 1) THEN
    G = GRAD (2, 1) **2 + GRAD (2, 2) **2
    IF ( N_GEOM == 3) G = G + GRAD (2, 3) **2
    F = GRAD (1, 1) * GRAD (2, 1) + GRAD (1, 2) * GRAD (2, 2)
    IF ( N_GEOM == 3) F = F + GRAD (1, 3) * GRAD (2, 3)
  END IF
  H = SQRT (E * G - F * F) ; ROOT_G = SQRT (G) ; ROOT_E = SQRT (E)
END SUBROUTINE SURFACE_LAME
                     
SUBROUTINE SURFACE_METRICS (P_GRAD, METRIC, E, F, G, H)
! * * * * * * * * * * * * * * * * * * * * * * * * * *
! COMPUTE SURFACE METRIC MEASURES & LAME PARAMETERS
!         METRIC = P_GRAD*P_GRAD_TRANPOSE
! * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants
 IMPLICIT NONE
 REAL(DP), INTENT(IN)  :: P_GRAD (N_PARM, N_SPACE)
 REAL(DP), INTENT(OUT) :: METRIC (N_PARM, N_PARM), E, F, G, H
                                                             
! N_SPACE = Dimension of physical space, >= N_PARM
! N_PARM  = Dimension of parametric space
! P_GRAD  = Parametric gradient of physical coordinates
! METRIC = Surface differential geometry metric measure
! METRIC = [ E  F ], Lame Parameters E, F, G, H
!          [ F  G ]  H**2 = E*G - F**2 = Det|Metric| > 0
                                                        
  METRIC = MATMUL (P_GRAD, TRANSPOSE(P_GRAD))
  E = METRIC (1, 1)
  F = METRIC (1, N_PARM)
  G = METRIC (N_PARM, N_PARM)
  H = SQRT (E * G - F**2)
END SUBROUTINE SURFACE_METRICS
                     
SUBROUTINE SURF_NORM_VEC (RJ, FFM, VECTOR)
! * * * * * * * * * * * * * * * * * * * * * * * * * *
! FIND UNIT NORMAL VECTOR TO A 3-D SURFACE AT A POINT
! * * * * * * * * * * * * * * * * * * * * * * * * * *
Use Precision_Module
 IMPLICIT NONE
 REAL(DP), INTENT(IN)  :: RJ (2, 3), FFM
 REAL(DP), INTENT(OUT) :: VECTOR (3)
 REAL(DP) :: ROOT
                                 
! VECTOR = GLOBAL COMPONENTS OF UNIT NORMAL VECTOR
                                                  
! THE FOLLOWING QUANTITIES ARE OBTAINED IN SUBR FUND_MAG_ONE
! RJ     = REDUCED JACOBIAN MATRIX
! A      = 1ST FUND MAGNITUDES(METRIC TENSOR), FFM = DET(A)
                                                           
  ROOT       = SQRT (FFM)
  VECTOR (1) = (RJ (1, 2)*RJ (2, 3) - RJ (1, 3)*RJ (2, 2))/ ROOT
  VECTOR (2) = (RJ (1, 3)*RJ (2, 1) - RJ (1, 1)*RJ (2, 3))/ ROOT
  VECTOR (3) = (RJ (1, 1)*RJ (2, 2) - RJ (1, 2)*RJ (2, 1))/ ROOT
END SUBROUTINE SURF_NORM_VEC
                            
SUBROUTINE TANGENT_VECTOR (DELTA, COORD, GRAD)
! * * * * * * * * * * * * * * * * * * * * * * * * *
!  CALCULATE THE TANGENT VECTORS AT A LOCAL POINT
!   ON A PARAMETRIC CURVE OR PARAMETRIC SURFACE
! * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants
 IMPLICIT NONE
 REAL(DP), INTENT(IN)  :: DELTA (N_PARM, N_GEOM)
 REAL(DP), INTENT(IN)  :: COORD (N_GEOM, N_SPACE)
 REAL(DP), INTENT(OUT) :: GRAD (N_PARM, N_SPACE)
                                                
! N_GEOM  = NUMBER OF NODES PER ELEMENT
! N_PARM  = NUMBER OF PARAMETRIC DIMENSIONS <= N_SPACE
! N_SPACE = DIMENSION OF PHYSICAL SPACE
! DELTA   = LOCAL  DERIVATIVES OF N_GEOM INTERPOLATION
!              FUNCTIONS AT POINT OF INTEREST.
! COORD   = SPATIAL COORDINATES OF ELEMENT'S NODES
! GRAD    = TANGENT MATRIX = DELTA*COORD
!              ROW 1 IS DR/DU, 2 DR/DV, 3 DR/DW
!              FOR N_PARM = N_SPACE, GRAD = JACOBIAN
    GRAD = MATMUL ( DELTA, COORD )
END SUBROUTINE TANGENT_VECTOR
                     
SUBROUTINE WHERE_IN_4Q (XY, COORD, INSIDE, R, S)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
! FIND LOCAL (R,S) COORDINATES OF POINT XY IN A Q4 ELEMENT
!   REF: C. HUA, FE IN ANAL & DESIGN, 7, 159-166, 1990
!            (SEE INSIDE_2D_POLYGON AS ALTERNATE)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants !b
!b Use Precision_Module
 IMPLICIT NONE
 REAL(DP), PARAMETER :: ZERO = 0.D0, HALF = 0.5D0, &
                        ONE = 1.D0, FOUR = 4.D0, TOL = 1.D-7
                                                             
 REAL(DP), INTENT(IN)  :: XY (2), COORD (4, 2)
 INTEGER,  INTENT(OUT) :: INSIDE
 REAL(DP), INTENT(OUT) :: R, S
  	
 REAL(DP) :: A, B, C, ROOT, R1, R2
 REAL(DP) :: TOABC (3, 4), ABC (3, 2), A1, A2, B1, B2, C1, C2
 REAL(DP) :: AB, AC, AD, BC, BD, DB, DC, D1, D2 ,TWO_A
 REAL(DP) :: ERROR_X, ERROR_Y
                             
! NODES ARE ASSUMED TO BE NUMBERED COUNTERCLOCKWISE
! COORD  = PHYSICAL COORDINATES OF NODES
! INSIDE = 0 IF XY IS NOT IN ELEMENT, ELSE = 1
! R,S    = NATURAL COORDINATES OF THE PT (-1 TO +1)
! XY     = PHSYICAL COORDINATES OF POINT INSIDE ELEMENT
!
  DATA TOABC / 1.D0, -1.D0, -1.D0,     -1.D0,  1.D0, -1.D0, &
               1.D0,  1.D0,  1.D0,     -1.D0, -1.D0,  1.D0 /
  EQUIVALENCE (ABC(1, 1), A1), (ABC(2, 1), B1), (ABC(3, 1), C1),&
              (ABC(1, 2), A2), (ABC(2, 2), B2), (ABC(3, 2), C2)
                                                               
!  INITALIZE AS OUTSIDE (ON EDGE IS INSIDE)
  INSIDE = 0  ; R = HUGE(R) ; S = HUGE(S)

!  LOOP OVER 4 SIDES TESTING TRIANGULAR AREA (+ IF IN)
  TWO_A = COORD (1, 1) * (COORD (2, 2) - XY (2) ) &
       + COORD (2, 1) * (XY (2) - COORD (1, 2) ) &
       + XY (1) * (COORD (1, 2) - COORD (2, 2) )
  IF ( TWO_A < -TOL) RETURN
    TWO_A = COORD (2, 1) * (COORD (3, 2) - XY (2) ) &
         + COORD (3, 1) * (XY (2) - COORD (2, 2) ) &
         + XY (1) * (COORD (2, 2) - COORD (3, 2) )
    IF ( TWO_A < -TOL) RETURN
      TWO_A = COORD (3, 1) * (COORD (4, 2) - XY (2) ) &
           + COORD (4, 1) * (XY (2) - COORD (3, 2) ) &
           + XY (1) * (COORD (3, 2) - COORD (4, 2) )
      IF ( TWO_A < -TOL) RETURN
        TWO_A = COORD (4, 1) * (COORD (1, 2) - XY (2) ) &
             + COORD (1, 1) * (XY (2) - COORD (4, 2) ) &
             + XY (1) * (COORD (4, 2) - COORD (1, 2) )
        IF ( TWO_A < -TOL) RETURN

!          POINT IS INSIDE, GET LOCAL COORDINATES
  INSIDE = 1
  D1 = 4. * XY (1) - COORD (1, 1) - COORD (2, 1) &
     - COORD (3, 1) - COORD (4, 1)
  D2 = 4. * XY (2) - COORD (1, 2) - COORD (2, 2) &
     - COORD (3, 2) - COORD (4, 2)
  ABC = MATMUL ( TOABC, COORD ) 
  AB = A1 * B2 - A2 * B1   ; AC = A1 * C2 - A2 * C1
  AD = A1 * D2 - A2 * D1   ; DC = D1 * C2 - D2 * C1

!          CHECK CASES
  IF ( (A1 * A2 * AB * AC) /= ZERO          .OR. &
       (A1 == ZERO .AND. (A2 * C1) /= ZERO) .OR. &
       (A2 == ZERO .AND. (A1 * B2) /= ZERO) ) THEN
    A = AB  ; B = C1 * B2 - C2 * B1 - AD  ; C = DC
    ROOT = SQRT (B * B - FOUR * A * C)
    R1 = HALF * (-B + ROOT) / A ; R2 = HALF * (-B - ROOT) / A
    R = R2
    IF ( -ONE <= R1 .AND. ONE >= R1 ) R = R1
    S = (AD-AB * R) / AC
  ELSEIF ( ( (A1 * A2) /= ZERO) .AND. (AB == ZERO) ) THEN
    R = A1 * DC / (B1 * AC + A1 * AD) ; S = AD / AC
  ELSEIF ( ( (A1 * A2) /= ZERO) .AND. (AC == ZERO) ) THEN
    DB = D1 * B2 - D2 * B1
    R = AD / AB  ; S = A1 * DB / (C1 * AB + A1 * AD)
  ELSE
    BC = B1 * C2 - B2 * C1   ; BD = B1 * D2 - B2 * D1
    R = DC / (A1 * D2 + BC)  ; S = BD / (A2 * D1 + BC)
  END IF

!               CHECK FOR ERRORS  
  ERROR_X = ABS (B1 * R + C1 * S - D1 + A1 * R * S)
  ERROR_Y = ABS (B2 * R + C2 * S - D2 + A2 * R * S)
  IF ( ERROR_X > TOL .OR. ERROR_Y > TOL ) WRITE (6, * ) &
      'WARNING: WHERE_IN_4Q ERROR_S ARE ', ERROR_X, ERROR_Y
  N_WARN = N_WARN + 1 ! INCREMENT WARNING
  IF ( -ONE > R .OR. ONE < R ) STOP 'WHERE_IN_4Q ERROR IN R'
  IF ( -ONE > S .OR. ONE < S ) STOP 'WHERE_IN_4Q ERROR IN S'
END SUBROUTINE WHERE_IN_4Q

SUBROUTINE TRI_GEOM (CODE, COORD, VALUE) 
! * * * * * * * * * * * * * * * * * * * * * * * * * *               
! EXACT GEOMETRIC PROPERTIES OF AN ARBITRARY TRIANGLE               
! * * * * * * * * * * * * * * * * * * * * * * * * * *               
Use Precision_Module
IMPLICIT NONE
  INTEGER,  INTENT (IN)  :: CODE
  REAL(DP), INTENT (IN)  :: COORD (3, 2) 
  REAL(DP), INTENT (OUT) :: VALUE
  REAL(DP), PARAMETER    :: PI = 3.1415926536d0
  REAL(DP) :: XI, XJ, XK, XB, YI, YJ, YK, YB, AREA

!  VALUE = INTEGRAL OF (X**M)*(Y**N)*DA (0<=M,N<=2)                  
!  VALUE BY CODE: 
!  1-AREA,                    2-VOL REV ABOUT Y AXIS,
!  3-VOL REV ABOUT Y,         4-FIRST MOMENT ABOUT X,
!  5-FIRST MOMENT ABOUT Y,    6-SEC MOMENT ABOUT X,                     
!  7-SEC MOMENT WRT X-Y AXES, 8-SEC MOMENT ABOUT Y                   

  XI = COORD (1, 1) ; XJ = COORD (2, 1) ; XK = COORD (3, 1) 
  YI = COORD (1, 2) ; YJ = COORD (2, 2) ; YK = COORD (3, 2) 
  AREA = 0.5 * (XI * (YJ - YK) + XJ * (YK - YI) + XK * (YI - YJ) ) 
  XB = (XI + XJ + XK) / 3.d0 ; YB = (YI + YJ + YK) / 3.d0 
      
  SELECT CASE (CODE)
    CASE (1) !     M=0, N=0,  AREA                             
      VALUE = AREA 
    CASE (2) !     M=0, N=0, VOLUME OF REVOLUTION ABOUT X-AXIS
      VALUE = 2.d0 * PI * AREA * YB 
      IF ( (ABS (YI) + ABS (YJ) + ABS (YK) ) /=             &
            ABS (YI + YJ + YK) ) WRITE (6, 5) ; 5 FORMAT    &
            ('TRI_GEOM: INVALID DATA FOR VOLUME OF REVOLUTION')
    CASE (3) !     M=0, N=0, VOLUME OF REVOLUTION ABOUT Y-AXIS
      VALUE = 2.d0 * PI * AREA * XB 
      IF ( (ABS (XI) + ABS (XJ) + ABS (XK) ) /=             &
            ABS (XI + XJ + XK) ) WRITE (6, 5) 
    CASE (4) !     M=0, N=1, FIRST MOMENT ABOUT X-AXIS
      VALUE = AREA * YB 
    CASE (5) !     M=1, N=0, FIRST MOMENT ABOUT Y-AXIS
      VALUE = AREA * XB 
    CASE (6) !     M=0, N=2, SECOND MOMENT ABOUT X-AXIS
      VALUE = AREA * (YI * YI + YJ * YJ + YK * YK &
            + 9.0 * YB * YB) / 12.d0
    CASE (7) !     M=1, N=1, SECOND MOMENT W.R.T. X-Y AXES
      VALUE = AREA * (XI * YI + XJ * YJ + XK * YK &
            + 9.0 * XB * YB) / 12.d0
    CASE (8) !     M=2, N=0, SECOND MOMENT ABOUT Y-AXIS
      VALUE = AREA * (XI * XI + XJ * XJ + XK * XK &
            + 9.0 * XB * XB) / 12.d0 
    CASE DEFAULT
      STOP 'TRI_GEOM: ARGUMENT OUT OF RANGE'
  END SELECT
END SUBROUTINE TRI_GEOM                         

SUBROUTINE TRI_CG_INERTIA (COORD, INERTIA) 
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * *         
! EXACT CENTROIDAL INERTIA TENSOR OF AN ARBITRARY TRIANGLE           
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * *         
Use Precision_Module
IMPLICIT NONE
integer, save :: count = 0
  REAL(DP), INTENT (IN)  :: COORD   (3, 2) 
  REAL(DP), INTENT (OUT) :: INERTIA (2, 2)
  REAL(DP) :: XI, XJ, XK, XB, YI, YJ, YK, YB, AREA

  XI = COORD (1, 1) ; XJ = COORD (2, 1) ; XK = COORD (3, 1) 
  YI = COORD (1, 2) ; YJ = COORD (2, 2) ; YK = COORD (3, 2) 
  AREA = 0.5 * (XI * (YJ - YK) + XJ * (YK - YI) + XK * (YI - YJ) ) 
  XB = (XI + XJ + XK) / 3.d0 ; YB = (YI + YJ + YK) / 3.d0 
      
!  With respect to global origin
  INERTIA (1, 1) =  AREA * (YI * YI + YJ * YJ + YK * YK &
                 + 9.0 * YB * YB) / 12.d0
  INERTIA (2, 1) = -AREA * (XI * YI + XJ * YJ + XK * YK &
                 + 9.0 * XB * YB) / 12.d0
  INERTIA (2, 2) =  AREA * (XI * XI + XJ * XJ + XK * XK &
                 + 9.0 * XB * XB) / 12.d0 
if ( count == 0 ) then
print *,'INERT ', AREA, XB, YB
print *,'INERT ', INERTIA
end if

!  With respect to centroid
  INERTIA (1, 1) = INERTIA (1, 1) - AREA * YB * YB
  INERTIA (2, 1) = INERTIA (2, 1) + AREA * XB * YB
  INERTIA (2, 2) = INERTIA (2, 2) - AREA * XB * XB
  INERTIA (1, 2) = INERTIA (2, 1)
if ( count == 0 ) then
count = 1                       
print *,'CG    ', INERTIA
end if
END SUBROUTINE TRI_CG_INERTIA                         

