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 3-dimensional
# triangular isoparametric membrane element with 3 nodes.
#
# -------------------------------------------------------------------

  persistent first = 1;
  persistent ng; persistent r; persistent  s; persistent w;
  persistent iu; persistent iv; persistent iw;  % Dof indices

# Initialize

  if (first)
     first = 0;
     ng = 3;               % cf. Bathe
     f  = 1 / 6;
     r  = [0.5, 0.5, 0];
     s  = [0, 0.5, 0.5];
     w  = f * [1, 1, 1];
     iu = 1 : 3 : 7; iv = iu + 1; iw = iv + 1;
  end

# Mass density

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

# Transformation to element coordinate system

  [x, y, TE] = mfs_esys(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

# Mass matrix in element coordinate system

  mE = sparse(9, 9);

  mE(iu, iu) = m33; mE(iv, iv) = m33; mE(iw, iw) = m33;

# Mass matrix in global coordinate system

  m = TE' * mE * TE;

end
