function BSE = mfs_BSE_s4(x, y)

# usage: BSE = mfs_BSE_s4(x, y)
#
# Input  x(4), y(4)  Arrays with nodal point coordinates
#
# Output BSE(4, 12)  Matrix to compute transverse shear strains
#                    at edge midpoints A, B, C and D from
#                    bending displacements w, phix, phiy
#
# The rows of matrix BSE correspond to points A, B, C and D.
#
# -------------------------------------------------------------------

  ind = [2, 3, 4, 1];

  BSE = zeros(4, 12);

# Differences of coordinates

  dx = 0.5 * (x(ind) - x); dy = 0.5 * (y(ind) - y);

# Point A (Edge 1-2)

  row = [-1, -dy(1), dx(1), 1, -dy(1), dx(1)];
  BSE(1, 1 : 6) = 0.25 * row;

# Point B (Edge 2-3)

  row = [-1, -dy(2), dx(2), 1, -dy(2), dx(2)];
  BSE(2, 4 : 9) = 0.25 * row;

# Point C (Edge 4-3)

  row = [1, dy(3), -dx(3), -1, dy(3), -dx(3)];
  BSE(3, 7 : 12) = 0.25 * row;

# Point D (Edge 1-4)

  row = [-1, dy(4), -dx(4), 1, dy(4), -dx(4)];
  BSE(4, [1 : 3, 10 : 12]) = 0.25 * row;

end
