function [BS, BK, JE] = mfs_BEE_s4(r, s, x, y)

# usage: [BS, BK, JE] = mfs_BEE_s4(r, s, x, y)
#
# Input  r, s     Parametric coordinates
#        x, y     Arrays with nodal point coordinates
#
# Output BS       Strain-displacement matrix
#        BK       Curvature-displacement matrix
#        JE       Jacobi determinant
#
# The function computes the strain-displacement matrix and the
# curvature-displacement matrix of a linear quadrilateral shell
# element with 4 nodes. The strain-displacement matrix does not
# contain the shear strain.
#
# -------------------------------------------------------------------

  persistent colx;  persistent coly;    % u, v
  persistent colrx; persistent colry;   % phi, psi

  if (isempty(colx))
     colx  = 1 : 2 : 7;  coly  = colx + 1;
     colrx = 2 : 3 : 11; colry = colrx + 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

  BS(1, colx) = dn(1, :); BS(2, coly) = dn(2, :);

# Curvature-displacement matrix

  BK([1, 3], colry) = dn; BK([2, 3], colrx) = -dn([2, 1], :);

end
