function [BM, BB, BS, JE] = mfs_BE3(r, s, x, y, reduc)

# usage: [BM, BB, BS, JE] = mfs_BE3(r, s, x, y, reduc)
#
# Input  r, s     Parametric coordinates
#        x, y     Arrays with nodal point coordinates
#        reduc    0 = return complete BM matrix
#                 1 = return reduced BM matrix (without shear)
#
# Output BM       Strain-displacement matrix for membrane strains
#        BB       Strain-displacement matrix for bending strains
#        BS       Matrix to compute transverse shear strain at
#                 integration point from transverse shear strain at
#                 edge midpoints
#        JE       Jacobi determinant
#
# The function computes the strain-displacement matrices of a linear
# quadrilateral shell element with 4 nodes.
#
# -------------------------------------------------------------------

  persistent first = 1;
  persistent cu;   persistent cv;     % u, v
  persistent cphi; persistent cpsi;   % phi, psi

  if (first)
     first = 0;
     cu   = 1 : 2 : 7; cv   = cu + 1;
     cphi = 1 : 2 : 7; cpsi = cphi + 1;  
  end

# Derivatives of shape functions with respect to r, s

  sm = 1 - s;  sp = 1 + s;
  rm = 1 - r;  rp = 1 + r;

  nr = 0.25 * [-sm, sm, sp, -sp];
  ns = 0.25 * [-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 for membrane strains

  if (reduc)
     BM(2, cv) = dn(2, :); BM(1, cu) = dn(1, :);
  else
     BM([3, 2], cv) = dn; BM([1, 3], cu) = dn;
  end

# Strain-displacement matrix for bending strains

  BB([1, 3], cpsi) = dn; BB([3, 2], cphi) = -dn;

# Strain-displacement matrix for transverse shear

  BS = Jinv * [sm, 0, sp, 0; 0, rp, 0, rm];
  
end
