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

# usage: [BE, JE] = mfs_BE_t6(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 quadratic
# triangular isoparametric element with 6 nodes.
#
# -------------------------------------------------------------------

  persistent colx; persistent coly;

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

# Derivatives of shape functions with respect to r, s

%  nr = zeros(1, 6); ns = zeros(1, 6);

  nr(4 : 6) = 4 * [1 - 2 * r - s, s, -s];
  ns(4 : 6) = 4 * [-r, r, 1 - r - 2 * s];

  nr(1 : 3) = [ -1 - 0.5 * (nr(4) + nr(6)),
                 1 - 0.5 * (nr(4) + nr(5)),
                   -0.5 * (nr(5) + nr(6))];
  ns(1 : 3) = [ -1 - 0.5 * (ns(4) + ns(6)),
                    -0.5 * (ns(4) + ns(5)),
                 1 - 0.5 * (ns(5) + ns(6))];

# 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(3, 12);
  BE([3, 2], coly) = dn; BE([1, 3], colx) = dn;

end
