function m = mfs_mle(elem)

# usage: m = mfs_mle(elem)
#
# Input  elem     Structure containing data of one element
#
# Output m        Element mass matrix
#
# The function computes the lumped 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 = 4;
     f  = 1 / sqrt(3);
     r  = [ -f, f, -f, f];
     s  = [ -f, -f, f, f];
     w  = [ 1, 1, 1, 1];
  end

# Mass density

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

# Area

  A    = 0;
  coor = elem.coor;

  for n = 1 : ng
      JE = mfs_JE_q9(r(n), s(n), coor(:, 1), coor(:, 2));
      A += w(n) * JE;
  end

# Mass matrix

  mp = rho * elem.geom.t * A / 16;
  m  = mp * speye(18); 
  m( 9 : 16,  9 : 16) = 2 * m( 9 : 16,  9 : 16);
  m(17 : 18, 17 : 18) = 4 * m(17 : 18, 17 : 18);

end
