function kE = mfs_ke(elem)

# usage: kE = mfs_ke(elem)
#
# Input  elem     Structure containing data of one element
#        id       Element identifier
#
# Output kE       Element stiffness matrix
#
# The function computes the stiffness matrix of a linear 3-dimensional
# triangular membrane element with 3 nodes with drilling dofs.
#
# ------------------------------------------------------------------------

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

  persistent ib;   % Indices of edge end nodes

# Initialize

  if (isempty(ng))
     ng = 3;
     f  = 1 / 6;
     r  = [0.5, 0.5, 0];
     s  = [0, 0.5, 0.5];
     w  = f * [1, 1, 1];
     ib = [2, 3, 1];
  end

# Material matrix and thickness

  E = mfs_emat2d(elem.mat);
  t = elem.geom.t;

# Transformation to element coordinate system

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

# Transformation matrix to t6 element

  T = mfs_rot2mid2d(coor) * TE;

# Add coordinates of midpoints

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

# Stiffness matrix of t6 element

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

  for n = 1 : ng
      [BE, JE] = mfs_BE_t6(r(n), s(n), coor(:, 1), coor(:, 2));
      wdV = w(n) * JE * t;
      kt6 = kt6 + BE' * (E * wdV) * BE;
  end

# Stiffness matrix of t6r element

  kE = T' * kt6 * T;

end
