function kE = mfs_ke(elem)

# usage: kE = mfs_ke(elem)
#
# Input  elem     Structure containing data of one element
#
# Output kE       Element stiffness matrix
#
# The function computes the stiffness matrix of a 2-dimensional 
# quadrilateral membrane element with 4 nodes and drilling dofs.
#
# -------------------------------------------------------------------

  persistent ng; persistent r;
  persistent  s; persistent w;

# Initialize Gauss integration

  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;
  end

# Material matrix, coordinates and thickness

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

# Transformation matrix to q4

  T = mfs_rot2mid2d(coor);

# Add coordinates of midpoints

  ib = [2, 3, 4, 1];
  coor = [coor; 0.5 * (coor + coor(ib, :))];

# Stiffness matrix of q8 element

  kq8 = zeros(16, 16,"double");

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

# Stiffness matrix of q8r element

  kE = T' * kq8 * T;

end
