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

# Mass of the element

  if (isfield(elem.mat, "rho"))
     d = elem.coor(2, :) - elem.coor(1, :);
     L = norm(d);
     rhoL = 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);

# Initialize mass matrix

  mE = zeros(12, 12, "double");

# Mass matrix for longitudinal vibrations

  ml  = (m0 / 6) * [2, 1; 1, 2];
  ixl = [ 1, 7];
  mE(ixl, ixl) = ml;

# Mass matrix for bending vibrations

  m11 = 156;
  m12 =  22 * L;
  m13 =  54;
  m14 = -13 * L;

  m22 =  4 * L^2;
  m23 = -m14;
  m24 = - 3 * L^2;

  m33 = 156;
  m34 = -m12;
  m44 = m22;

  mb = m0 / 420;

# Bending vibrations in xy-plane

  mby = mb  * [ m11, m12, m13, m14;
                m12, m22, m23, m24;
                m13, m23, m33, m34;
                m14, m24, m34, m44 ];

  ixby = [ 2, 6, 8, 12];
  mE(ixby, ixby) = mby;

# Bending vibrations in xz-plane

  mbz = mb  * [ m11, -m12, m13, -m14;
                -m12, m22, -m23, m24;
                m13, -m23, m33, -m34;
                -m14, m24, -m34, m44 ];

  ixbz = [ 3, 5, 9, 11];
  mE(ixbz, ixbz) = mbz;

# Mass matrix for torsional vibrations

  mt  = (Jx / 6) * [2, 1; 1, 2];
  ixt = [4, 10];
  mE(ixt, ixt) = mt;

# 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' * mE * T;
  
end
