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

# usage: [BE, JE] = mfs_BEE_q4(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 linear
# quadrilateral isoparametric element with 4 nodes. The matrix does
# not include the shear strains.
#
# -------------------------------------------------------------------

  persistent colx; persistent coly;

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

# Derivatives of shape functions with respect to r, s

  sm = 0.25 * (1 - s);  sp = 0.25 * (1 + s);
  rm = 0.25 * (1 - r);  rp = 0.25 * (1 + r);

  nr = [-sm, sm, sp, -sp];
  ns = [-rm, -rp, rp, rm];

# 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(2, 8);
  BE(1, colx) = dn(1, :); BE(2, coly) = dn(2, :);

end
