function JE = mfs_JE_q8(r, s, x, y)

# usage: JE = mfs_JE_q8(r, s, x, y)
#
# Input  r, s     Parametric coordinates
#        x, y     Arrays with nodal point coordinates
#
# Output JE       Jacobi determinant
#
# The function computes the Jacobi determinant of a quadratic
# quadrilateral isoparametric element with 8 nodes.
#
# -------------------------------------------------------------------

# Derivatives of shape functions with respect to r, s

%  nr = zeros(1, 8); ns = zeros(1, 8);

  sm = 1 - s;  sp = 1 + s;  rm = 1 - r;  rp = 1 + r;
  sq = sm * sp;  rq = rm * rp;

  nr(5 : 8) = [ -r * sm, 0.5 * sq, -r * sp, -0.5 * sq];
  ns(5 : 8) = [ -0.5 * rq, -s * rp, 0.5 * rq, -s * rm];

  sm = 0.25 * sm;  sp = 0.25 * sp;
  rm = 0.25 * rm;  rp = 0.25 * rp;

  nr(1 : 4) = [ -sm - 0.5 * (nr(5) + nr(8)),
                 sm - 0.5 * (nr(5) + nr(6)),
                 sp - 0.5 * (nr(6) + nr(7)),
                -sp - 0.5 * (nr(7) + nr(8))];
  ns(1 : 4) = [ -rm - 0.5 * (ns(5) + ns(8)),
                -rp - 0.5 * (ns(5) + ns(6)),
                 rp - 0.5 * (ns(6) + ns(7)),
                 rm - 0.5 * (ns(7) + ns(8))];

# Jacobi matrix

  Jmat = [nr; ns] * [x, y];

# Jacobi determinant

  JE = abs(det(Jmat));

end
