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

  persistent first = 1;
  persistent ng; persistent r;
  persistent  s; persistent w;

# Initialize Gauss integration

  if (first)
     first = 0;
     ng = 4;
     f  = 1 / sqrt(3);
     r  = [ -f, f, -f, f];
     s  = [ -f, -f, f, f];
     w  = [ 1, 1, 1, 1];
  end

# Material matrix and thickness

  E   = mfs_emat2d(elem.mat);
  E22 = E([1 : 2], [1 : 2]);
  G   = E(3, 3);
  t   = elem.geom.t;

# Transformation to element coordinate system

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

# Stiffness matrix in element coordinate system

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

  for n = 1 : ng
      [BE, JE] = mfs_BEE_q4(r(n), s(n), x, y);
      wdV = w(n) * JE * t;
      kEE = kEE + BE' * (E22 * wdV) * BE;
  end

  [BEG, JE] = mfs_BEG_q4(0, 0, x, y);
  wdV = 4 * JE * t;
  kEE = kEE + BEG' * (G * wdV) * BEG;

# Stiffness matrix in global coordinate system

  kE = TE' * kEE * TE;

end
