function [BE, JE] = mfs_BE_q9(r, s, x, y)

# usage: [BE, JE] = mfs_BE_q9(r, s, x, y)
#
# Input  r, s     Parametric coordinates
#        x, y     Arrays with nodal point coordinates
#
# Output BE       Strain-displacement matrix
#        JE       Jacobi determinant
#
# The function computes the strain-displacement matrix of a quadratic
# quadrilateral isoparametric element with 9 nodes.
#
# -------------------------------------------------------------------

  persistent colx; persistent coly;

  if (isempty(colx))
     colx = 1 : 2 : 17; coly = colx + 1;
  end

# Derivatives of shape functions with respect to r, s

%  nr = zeros(1, 9); ns = zeros(1, 9);  % more efficient without

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

  nr9 = - r * sq; nr(9) = 2 * nr9;
  ns9 = - s * rq; ns(9) = 2 * ns9;

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

  sm  = 0.25 * sm;  sp = 0.25 * sp;
  rm  = 0.25 * rm;  rp = 0.25 * rp;
  nr9 = 0.5 * nr9;  ns9 = 0.5 * ns9;

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

# Jacobi matrix

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

# Inverse Jacobi matrix

  JE   = abs(det(Jmat));
  Jinv = [Jmat(2,2), - Jmat(1, 2);
          -Jmat(2, 1), Jmat(1, 1)] / JE;

# Derivatives of shape functions with respect to x, y

  dn = Jinv * nrs;

# Strain-displacement matrix

%  BE = zeros(3, 18);
  BE([3, 2], coly) = dn; BE([1, 3], colx) = dn;

end
