function [BM, BB, BS0, BS1, JE] = mfs_BE(x, y)

# usage: [BM, BB, BS0, BS1, JE] = mfs_BEMB(x, y)
#
# Input  x, y     Arrays with nodal point coordinates
#
# Output BM     Strain-displacement matrix for membrane strains
#        BB     Strain-displacement matrix for bending strains
#        BS0    Strain-displacement matrix for shear strains (const.)
#        BS1    Strain-displacement matrix for shear strains (linear)
#        JE     Jacobi determinant
#
# The function computes the strain-displacement matrices of a linear
# triangular shell element with 3 nodes.
#
# -------------------------------------------------------------------

  persistent first = 1;
  persistent mu;   persistent mv;    % u, v membrane
  persistent bphi; persistent bpsi;  % phi, psi bending
  persistent sw;                     % w shear
  persistent sphi; persistent spsi;  % phi, psi shear

  if (first)
     first = 0;
     mu   = 1 : 2 : 5; mv   = mu + 1;
     bphi = 1 : 2 : 5; bpsi = bphi + 1;
     sw   = 1 : 3 : 7; sphi = sw + 1; spsi = sphi + 1;
  end

# Differences of coordinates

  dx12 = x(2) - x(1); dx13 = x(3) - x(1); dx23 = x(3) - x(2);
  dy12 = y(2) - y(1); dy13 = y(3) - y(1); dy23 = y(3) - y(2);

# Jacobi determinant

  JE = dx12 * dy13 - dx13 * dy12;

# Inverse Jacobi matrix

  Jinv = [dy13, -dy12; -dx13, dx12] / JE;

# Shape function derivatives

  dn  = [-dy23,  dy13, -dy12;          % with respect to x
          dx23, -dx13,  dx12] / JE;    % with respect to y

# Strain-displacement matrix for membrane strains

  BM([3, 2], mv) = dn; BM([1, 3], mu) = dn;

# Strain-displacement matrix for bending strains

  BB([1, 3], bpsi) = dn; BB([3, 2], bphi) = -dn;

# Matrices for transverse shear strain interpolation

  A = 0.25 * [-2, -dy12, dx12, 2, -dy12, dx12, 0, 0, 0];
  B = 0.25 * [-2, -dy13, dx13, 0, 0, 0, 2, -dy13, dx13];
  C = 0.25 * [0, -dy23, dx23, 0, dy13, -dx13, 0, -dy12, dx12];

  BS0 = Jinv * [A; B];
  BS1 = Jinv * [C; -C];

end
