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 shell element with 4
# nodes.
#
# ------------------------------------------------------------------------

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

  persistent kS; persistent a0;

  persistent mdofs; persistent bdofs; persistent sdofs; persistent ddofs;

# Initialize

  if (first)
     first = 0;
     ng = 4;                % Gauss integration
     f  = 1 / sqrt(3);
     r  = [ -f, f, -f, f];
     s  = [ -f, -f, f, f];
     w  = [ 1, 1, 1, 1];
     kS = 5 / 6;            % Shear correction factor
     a0 = 1e-2;             % Artificial stiffness factor
     mdofs = [1, 2, 7, 8, 13, 14, 19, 20];                 % Membrane
     bdofs = [4, 5, 10, 11, 16, 17, 22, 23];               % Bending
     sdofs = [3, 4, 5, 9, 10, 11, 15, 16, 17, 21, 22, 23]; % Shear
     ddofs = [6, 12, 18, 24];                              % Drilling
  end

  if (isfield(elem.geom, "ka"))
     a = elem.geom.ka;
  else
     a = a0;
  end

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

# Material matrix and thickness

  h   = elem.geom.t;
  E   = mfs_emat2d(elem.mat);
  E_m = E(1 : 2, 1 : 2) * h;
  G   = E(3, 3) * h;
  Gs  = kS * G;
  E_b = (h^3/12) * E;

# Transformation to element coordinate system

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

# Initialize bending and shear contributions

  k_b = zeros(8, 8, "double");
  k_s = zeros(12, 12, "double");

# Contribution of shear to membrane part

  [BM, JE] = mfs_BEG_q4(0, 0, x, y);
  k_m      = 4 * JE * G * BM' * BM;

# Ensure convergence to Kirchhoff theory

  Gs = mfs_Gscale(Gs, 4 * JE, elem.geom);

# Transverse shear at edge points

  BSE = mfs_BSE_s4(x, y);

# Integration point loop for membrane, bending and transverse shear

  for n = 1 : ng
      [BM, BB, BG, JE] = mfs_BE3(r(n), s(n), x, y, 1);
      BS = BG * BSE;
      wdA = w(n) * JE;
      k_m = k_m + BM' * (E_m * wdA) * BM;
      k_b = k_b + BB' * (E_b * wdA) * BB;
      k_s = k_s + (Gs  * wdA) * BS' * BS;
  end

# Drilling dof stabilization

  kbm = mean(diag(k_b));
  b   = a * kbm; bd = 2 * b;
  k_d = [bd, -b, 0, -b;
         -b, bd, -b, 0;
         0, -b, bd, -b;
         -b, 0, -b, bd];

# Complete stiffness matrix in element coordinate systems

  kEE(mdofs, mdofs)  = k_m;
  kEE(bdofs, bdofs)  = k_b;
  kEE(sdofs, sdofs) += k_s;
  kEE(ddofs, ddofs)  = k_d;

# Stiffness matrix in global coordinate system

  kE = TE' * kEE * TE;

end
