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
# quadrilateral isoparametric membrane element with 9 nodes.
#
# -------------------------------------------------------------------

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

# Initialize Gauss integration

  if (isempty(ng))
     ng = 9;
     f = sqrt(3/5);
     r = [ -f, 0, f, -f, 0, f, -f, 0, f];
     s = [ -f, -f, -f, 0, 0, 0, f, f, f];
     w = [ 25, 40, 25, 40, 64, 40, 25, 40, 25] / 81;
  end

# Initialize result

  m = sparse(18, 18);

# Mass density

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

# Degree of freedom indices

  u = 1 : 2 : 17;
  v = 2 : 2 : 18; 

# Integration loop

  m99 = zeros(9, 9, "double");

  coor = elem.coor;
  for n = 1 : ng
      JE = mfs_JE_q9(r(n), s(n), coor(:, 1), coor(:, 2));
      wdA = w(n) * JE;
      N   = mfs_NE_q9(r(n), s(n));
      m99 += N' * N * wdA;
  end

  m99     = rho * elem.geom.t * m99;
  m(u, u) = m99;
  m(v, v) = m99;

end
