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 3
# nodes.
#
# ------------------------------------------------------------------------

  persistent first = 1;
  persistent ng; persistent r; persistent  s; persistent w;
  persistent kS; persistent a0;
  persistent sw;                     % w shear
  persistent sphi; persistent spsi;  % phi, psi shear
  persistent mdofs; persistent bdofs; persistent sdofs; persistent ddofs;

# Initialize

  if (first)
     first = 0;
     ng = 3;               % cf. Bathe
     w  = 1 / 6;
     r  = [0.5, 0.5, 0];
     s  = [0, 0.5, 0.5];
     kS = 5 / 6;            % Shear correction factor
     a0 = 1e-2;             % Artificial stiffness factor
     sw = 1 : 3 : 7; sphi = sw + 1; spsi = sphi + 1;
     mdofs = [1, 2, 7, 8, 13, 14];             % Membrane
     bdofs = [4, 5, 10, 11, 16, 17];           % Bending
     sdofs = [3, 4, 5, 9, 10, 11, 15, 16, 17]; % Shear
     ddofs = [6, 12, 18];                      % Drilling
  end

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

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

# Material matrix and thickness

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

# Transformation to element coordinate system

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

# Strain-displacement matrices

  [BM, BB, BS0, BS1, JE] = mfs_BE(x, y);
  Area = 0.5 * JE;

# Membrane part of stiffness matrix

  k_m = BM' * (Area * E_m) * BM;

# Bending part of stiffness matrix

  k_b = BB' * (Area * E_b) * BB;

# Transverse shear part of stiffness matrix

  k_s = zeros(9, 9);
  Gs  = mfs_Gscale(Gs, JE, elem.geom);
  AG  = w * JE * Gs;

  for n = 1 : ng
      BS  = 2 * (BS0 + [s(n); r(n)] .* BS1);
      k_s += AG * BS' * BS;
  end

# Drilling dof stabilization

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

# Complete stiffness matrix in element coordinate system

  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
