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 3-dimensional beam 
# element.
#
# -------------------------------------------------------------------

# Mass of the element

  if (isfield(elem.mat, "rho"))
     d = elem.coor(2, :) - elem.coor(1, :);
     L = norm(d);
     rhoL = 0.5 * L * elem.mat.rho;
     m0   = rhoL* elem.geom.A;
     Jx   = rhoL * (elem.geom.Iy + elem.geom.Iz);
  else
     m = sparse(12, 12);
     return;
  end

# Element coordinate system

  if (isfield(elem.geom, "Q"))
     v = elem.geom.Q - elem.coor(1, :);
  else
     v = elem.geom.v;
  end

  ex = d' / L;
  ez = v' - (v * ex) * ex;
  Lz = norm(ez);
  ez = ez / Lz;
  ey = cross(ez, ex);

# Mass matrix in element coordinate system

  m = diag([m0, m0, m0, Jx, 0, 0, m0, m0, m0, Jx, 0, 0]);

# Transformation matrix from structural to element system

  T33 = [ex'; ey'; ez'];
  T   = zeros(12, 12, "double");
  T( 1 :  3,  1 : 3) = T33;
  T( 4 :  6,  4 : 6) = T33;
  T( 7 :  9,  7 : 9) = T33;
  T(10 : 12, 10 :12) = T33;

# Transformation matrix to center of area

  if (isfield(elem.geom, "P"))
     TA = eye(12, 12);
     TA(1, 5) = -elem.geom.P(2);
     TA(1, 6) =  elem.geom.P(1);
     TA(2, 4) = -TA(1, 5);
     TA(3, 4) = -TA(1, 6);
     TA(7, 11) = -elem.geom.P(2);
     TA(7, 12) =  elem.geom.P(1);
     TA(8, 10) = -TA(7, 11);
     TA(9, 10) = -TA(7, 12);
     T = TA * T;
  end

# Mass matrix in structure coordinate system

  m = T' * m * T;
  
end
