getJacobian_2D Module Function

  module function getJacobian_2D(N, xi, eta, xy, alpha) result(J)
      !*
      ! Calculates the Jacobian of a quadrilateral element at a point
      !
      ! The Jacobian of an element is defined as:
      ! \[ \boldsymbol{J} = \boldsymbol{P} \boldsymbol{X} \]
      !

      ! Where:
      ! \[ \boldsymbol{P} = \left[ \begin{array}{cc}
      !     \frac{\partial H_1}{\partial \xi} & \frac{\partial H_2}{\partial \xi} \\
      !     \frac{\partial H_1}{\partial \eta} & \frac{\partial H_2}{\partial \eta} \end{array}
      !       \cdots
      !     \begin{array}{cc}
      !     \frac{\partial H_{N-1}}{\partial \xi} & \frac{\partial H_{N}}{\partial \xi} \\
      !     \frac{\partial H_{N-1}}{\partial \eta} & \frac{\partial H_{N}}{\partial \eta} \end{array}
      ! \right]\]
      !
      ! \[ \boldsymbol{X} = \left[ \begin{array}{c}
      !      \begin{array}{cc}
      !        x_1 & y_1 \\
      !        x_2 & y_2
      !      \end{array} \\\\
      !      \vdots \\\\
      !      \begin{array}{cc}
      !        x_{N-1} & y_{N-1} \\
      !        x_N & y_N
      !      \end{array} \\
      !    \end{array} \right]\]
      !
      ! Thus:
      !
      ! \[ \boldsymbol{J} =  \left[ \begin{array}{cc}
      !           \boldsymbol{H}_{\xi} \cdot \boldsymbol{x} & \boldsymbol{H}_{\xi} \cdot \boldsymbol{y} \\
      !           \boldsymbol{H}_{\eta} \cdot \boldsymbol{x} & \boldsymbol{H}_{\eta} \cdot \boldsymbol{y} \\
      !         \end{array} \right] = \left[ \begin{array}{cc}
      !           \frac{\partial \xi}{\partial x} & \frac{\partial \xi}{\partial y} \\
      !           \frac{\partial \eta}{\partial x} & \frac{\partial \eta}{\partial y} \\
      !         \end{array} \right] \]
      !
      ! * \( \boldsymbol{H} \) : Vector of basis functions
      ! * \( \boldsymbol{x} \) : Vector of X-coordinates for all element nodes
      ! * \( \boldsymbol{y} \) : Vector of Y-coordinates for all element nodes
      !
      ! * \( H_i \) : Basis function \(i\)
      ! * \( x_i \) : X-coordinate of node \(i\)
      ! * \( y_i \) : Y-coordinate of node \(i\)

      integer,                  intent(in)  :: N      !! Number of points in element
      real(wp),                 intent(in)  :: xi     !! Value of xi coordinate in isoparametric element
      real(wp),                 intent(in)  :: eta    !! Value of eta coordinate in isoparametric element
      real(wp), dimension(N,2), intent(in)  :: xy     !! Array of x and y coordinates for the element
      real(wp), dimension(N,N), intent(in)  :: alpha  !! Array of coefficients used to define basis functions
      real(wp), dimension(2,2)              :: J      !! Jacobian of the isoparametric element at (xi,eta)

      real(wp), parameter                   :: eps = epsilon(0e0)
      real(wp), dimension(2,N)              :: P      !! Array of derivatives of each basis function at (xi,eta)
      real(wp), dimension(N)                :: x      !! Row in element coefficient matrix

      ! P is a matrix containing derivatives of each basis function at (xi,eta)
      ! P = [dN_1/dxi, dN_2/dxi, dN_3/dxi, ...
      !      dN_1/deta, dN_2/deta, dN_3/deta, ...]

      P = 0._wp
      P(1,:) = [getArow(N, xi+eps, eta    ) - &
          getArow(N, xi-eps, eta    )]
      P(2,:) = [getArow(N, xi    , eta+eps) - &
          getArow(N, xi    , eta-eps)]

      P = P / ( 2._wp*eps )

      P = matmul(P,alpha)

      J = matmul(P,xy)

      return
  end function getJacobian_2D