! Common subroutines and parameters for mesh elements
module moduleMeshCommon
  ! values for integration gauss integral
  ! Segment
  real(8), parameter:: corSeg(1:3)  = (/ -dsqrt(3.D0/5.D0),  0.D0,      dsqrt(3.D0/5.D0) /)
  real(8), parameter:: wSeg(1:3)    = (/        5.D0/9.D0,   8.D0/9.D0,       5.D0/9.D0  /)

  ! Quad
  real(8), parameter:: corQuad(1:3) = (/ -dsqrt(3.D0/5.D0),  0.D0,      dsqrt(3.D0/5.D0) /)
  real(8), parameter:: wQuad(1:3)   = (/        5.D0/9.D0,   8.D0/9.D0,       5.D0/9.D0  /)

  ! tria
  real(8), parameter:: Xi1Tria(1:4) = (/        1.D0/3.D0,   1.D0/5.D0,       3.D0/5.D0,   1.D0/5.D0  /)
  real(8), parameter:: Xi2Tria(1:4) = (/        1.D0/3.D0,   1.D0/5.D0,       1.D0/5.D0,   3.D0/5.D0  /)
  real(8), parameter:: wTria(1:4)   = (/      -27.D0/96.D0, 25.D0/96.D0,     25.D0/96.D0, 25.D0/96.D0 /)

  contains
    ! ELEMENT FUNCTIONS
    ! Point
    pure function fPsiPoint(Xi, nNodes) RESULT(fPsi)
      implicit none

      real(8), intent(in):: Xi(1:3)
      integer, intent(in):: nNodes
      real(8):: fPsi(1:nNodes)

      fPsi = (/ 1.D0 /)

    end function fPsiPoint

    ! Segment
    pure function fPsiSegm(Xi, nNodes) result(fPsi)
      implicit none

      real(8), intent(in):: Xi(1:3)
      integer, intent(in):: nNodes
      real(8):: fPsi(1:nNodes)

      fPsi = (/ 1.D0 - Xi(1), &
                1.D0 + Xi(1) /)

      fPsi    = fPsi * 0.50D0

    end function fPsiSegm

    ! Quad
    pure function fPsiQuad(Xi, nNodes) result(fPsi)
      implicit none

      real(8), intent(in):: Xi(1:3)
      integer, intent(in):: nNodes
      real(8):: fPsi(1:nNodes)

      fPsi = 0.D0

      fPsi = (/ (1.D0 - Xi(1)) * (1.D0 - Xi(2)), &
                (1.D0 + Xi(1)) * (1.D0 - Xi(2)), &
                (1.D0 + Xi(1)) * (1.D0 + Xi(2)), &
                (1.D0 - Xi(1)) * (1.D0 + Xi(2)) /)

      fPsi    = fPsi * 0.25D0

    end function fPsiQuad

    ! Tria
    pure function fPsiTria(Xi, nNodes) result(fPsi)
      IMPLICIT NONE

      real(8), intent(in):: Xi(1:3)
      integer, intent(in):: nNodes
      real(8):: fPsi(1:nNodes)

      fPsi(1) = 1.D0 - Xi(1) - Xi(2)
      fPsi(2) =        Xi(1)
      fPsi(3) =                Xi(2)

    end function fPsiTria

    ! Tetra
    pure function fPsiTetra(Xi, nNodes) result(fPsi)
      implicit none

      real(8), intent(in):: Xi(1:3)
      integer, intent(in):: nNodes
      real(8):: fPsi(1:nNodes)

      fPsi(1) = 1.D0 - Xi(1) - Xi(2) - Xi(3)
      fPsi(2) =        Xi(1)
      fPsi(3) =                Xi(2)
      fPsi(4) =                        Xi(3)

    end function fPsiTetra

    ! DERIVATIVE FOR ELEMENT FUNCTIONS
    ! Segment
    pure function dPsiSegm(Xi, nNodes) result(dPsi)
      implicit none

      real(8), intent(in):: Xi(1:3)
      integer, intent(in):: nNodes
      real(8):: dPsi(1:3,1:nNodes)

      dPsi = 0.D0

      dPsi(1, 1:2) = (/ -5.D-1, 5.D-1 /)

    end function dPsiSegm

    ! Quad
    pure function dPsiQuad(Xi, nNodes) result(dPsi)
      implicit none

      real(8), intent(in):: Xi(1:3)
      integer, intent(in):: nNodes
      real(8):: dPsi(1:3,1:nNodes)

      dPsi = 0.D0

      dPsi(1, 1:4) = (/ -(1.D0 - Xi(2)), &
                         (1.D0 - Xi(2)), &
                         (1.D0 + Xi(2)), &
                        -(1.D0 + Xi(2)) /)

      dPsi(2, 1:4) = (/ -(1.D0 - Xi(1)), &
                        -(1.D0 + Xi(1)), &
                         (1.D0 + Xi(1)), &
                         (1.D0 - Xi(1)) /)

      dPsi = dPsi * 0.25D0

    end function dPsiQuad

    ! Tria
    pure function dPsiTria(Xi, nNodes) result(dPsi)
      implicit none

      real(8), intent(in):: Xi(1:3)
      integer, intent(in):: nNodes
      real(8):: dPsi(1:3,1:nNodes)

      dPsi = 0.D0

      dPsi(1,1:3) = (/ -1.D0, 1.D0, 0.D0 /)
      dPsi(2,1:3) = (/ -1.D0, 0.D0, 1.D0 /)

    end function dPsiTria

    !Compute element derivative functions in point Xi
    pure function dPsiTetra(Xi, nNodes) result(dPsi)
      implicit none

      real(8), intent(in):: Xi(1:3)
      integer, intent(in):: nNodes
      real(8):: dPsi(1:3, 1:nNodes)

      dPsi = 0.D0

      dPsi(1,1:4) = (/ -1.D0, 1.D0, 0.D0, 0.D0 /)
      dPsi(2,1:4) = (/ -1.D0, 0.D0, 1.D0, 0.D0 /)
      dPsi(3,1:4) = (/ -1.D0, 0.D0, 0.D0, 1.D0 /)

    end function dPsiTetra

end module moduleMeshCommon