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 4 nodes.
#
# -------------------------------------------------------------------

  persistent first = 1;
  persistent ib;    % Indices of edge end nodes
  persistent mdofs; % Membrane degrees of freedom

# Initiliaze

  if (first)
     first = 0;
     ib    = [2, 3, 4, 1];
     mdofs = [1 : 3, 7 : 9, 13 : 15, 19 : 21];
  end

  m = sparse(24, 24);

# Mass density

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

# Area

  coor  = elem.coor;
  edges = coor(ib, :) - coor;
  n     = cross(edges([1, 3], :), edges([2, 4], :));
  A     = mean(norm(n, "rows"));

# Mass matrix

  m(mdofs, mdofs) = 0.25 * rho * elem.geom.t * A * speye(12, 12);

end
