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 membrane element with 3 nodes and drilling dofs.
#
# -------------------------------------------------------------------

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

# Initialize Gauss integration (cf. Bathe)

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

# Mass density

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

# Transformation matrix to t6

  coor = elem.coor;
  T    = mfs_rot2mid2d(coor);

# Add coordinates of midpoints

  ib   = [2, 3, 1];
  coor = [coor; 0.5 * (coor + coor(ib, :))];

# Mass matrix of t6 element

  m66 = zeros(6, 6, "double");

  for n = 1 : ng
      JE = mfs_JE_t6(r(n), s(n), coor(:, 1), coor(:, 2));
      wdA = w(n) * JE;
      N   = mfs_NE_t6(r(n), s(n));
      m66 += N' * N * wdA;
  end

  m66 = rho * elem.geom.t * m66;

  u = 1 : 2 : 11; v = 2 : 2 : 12; 

  mt6 = zeros(12, 12, "double");
  mt6(u, u) = m66; mt6(v, v) = m66;

# Mass matrix of t6r element

  m = T' * mt6 * T;

end
