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

  persistent ng; % Gauss integration
  persistent r; persistent  s; persistent w;

  persistent ib; % Indices of edge end nodes
  persistent iu; persistent iv; persistent iw; % Dof indices
  persistent iuvr;
  persistent T84;

# Initialize

  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;
     ib = [2, 3, 4, 1];
     iu = [1 : 2 : 15]; iv = iu + 1;                   % refers to q8
     iw = [3, 7, 11, 15];                              % refers to q4r
     iuvr = [1, 2, 4, 5, 6, 8, 9, 10, 12, 13, 14, 16]; % refers to q4r
     T84 = [speye(4, 4); 
            0.5, 0.5, 0, 0; 0, 0.5, 0.5, 0;
            0, 0, 0.5, 0.5; 0.5, 0, 0, 0.5];  % q8 -> q4
  end

# Mass density

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

# Transformation to element coordinate system

  [coor, TE] = mfs_esys(elem.coor, 2, elem.geom);

# Transformation to q8 element

  TR =  mfs_rot2mid2d(coor);

# Add coordinates of midpoints

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

# Mass matrix of q8 element

  m88 = zeros(8, 8, "double");

  for n = 1 : ng
      JE = mfs_JE_q8(r(n), s(n), coor(:, 1), coor(:, 2));
      wdA = w(n) * JE;
      N   = mfs_NE_q8(r(n), s(n));
      m88 += N' * N * wdA;
  end

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

  mq8         = sparse(16, 16);
  mq8(iu, iu) = m88; mq8(iv, iv) = m88;

# Mass matrix of 2-dim. q4r element

  m2dq4r = TR' * mq8 * TR;

# Mass matrix of 3-dim. q4r element in local coordinate system

  ml  = zeros(16, 16, "double"); 
  m44 = T84' * m88 * T84;

  ml(iuvr, iuvr) = m2dq4r;
  ml(iw, iw)     = m44;

# Transformation to global coordinate system

  m = TE' * ml * TE;

end
