function m = mfs_mce(elem)

# usage: m = mfs_mce(elem)
#
# Input  elem     Structure containing data of one element
#
# Output m        Element mass matrix
#
# The function computes the consistent mass matrix of a 2-dimensional
# triangular isoparametric membrane element with 3 nodes.
#
# -------------------------------------------------------------------

  persistent ng; persistent r;
  persistent  s; persistent w;

# Initialize Gauss integration (cf. Bathe)

  if (isempty(ng))
     ng = 3;
     f  = 1 / 6;
     r  = f * [1, 4, 1];
     s  = f * [1, 1, 4];
     w  = 2 * f;
  end

# Initialize result

  m   = sparse(6, 6);

# Mass density

  if (isfield(elem.mat, "rho"))
     rho = elem.mat.rho;
  else
     return;
  end

# Degree of freedom indices

  u = 1 : 2 : 5;
  v = 2 : 2 : 6; 

# Coordinates

  x  = elem.coor(:, 1);  y = elem.coor(:, 2);

# Area and weighted mass

  A  = 0.5 * ((x(2) - x(1)) * (y(3) - y(1)) - (x(3) - x(1)) * (y(2) - y(1)));
  wm = rho * elem.geom.t * A / 3;

# Integration Loop

  m33 = zeros(3, 3, "double");

  for n = 1 : ng
      N   = [1 - r(n) - s(n), r(n), s(n)];
      m33 += wm * N' * N;
  end

  m(u, u) = m33;
  m(v, v) = m33;

end
