function [coore, TE, T33] = mfs_esys(coor, geom)

# usage: [coore, TE, T33] = mfs_esys(coor, geom)
#
# Input  coor(4, 3)   Coordinates of element nodes
#        geom         Geometrical data
#
# Output coore(4, 2)  Coordinates of element nodes in 2d element
#                     coordinate system
#        TE(24, 24)   Matrix to transform displacements and rotations
#                     to element coordinate system
#        T33(3, 3)    Elementary transformation matrix to element
#                     coordinate system
#
# The function computes the x- and y-coordinates of the element
# nodal points with respect to the element coordinate system and the
# transformation matrix from the global to the element coordinate
# system for quadrilateral shell elements.
#
# --------------------------------------------------------------------

# Origin

  C = mean(coor);                  % Centre of element

# Basis vectors of element coordinate system

  M23 = mean(coor([2, 3], :));     % Midpoint between points 2 and 3
  M34 = mean(coor([3, 4], :));     % Midpoint between points 3 and 4

  bx = M23 - C;
  bx = bx / norm(bx);
  c2 = M34 - C;
  bz = cross(bx, c2);
  bz = bz / norm(bz);
  by = cross(bz, bx);

  T33  = [bx; by; bz];

# Coordinates in element coordinate system

  for n = 1 : 4
      coor(n, :) = coor(n, :) - C;
  end

  coore = coor * [bx', by'];

# Transformation matrix

  TE = sparse(24, 24);

  T66 = sparse(6, 6);
  T66(1 : 3, 1 : 3) = T33; T66(4 : 6, 4 : 6) = T33;

  if (isfield(geom, "zo"))
     zo = geom.zo;
     TP = speye(6, 6);
     TP(1, 5) = -zo; TP(2, 4) = zo;
     T66 = TP * T66;
  end

  ix = 1 : 6;
  for n = 1 : 4
      TE(ix, ix) = T66;
      ix += 6; 
  end

# Check for warping

  x = coore(:, 1); y = coore(:, 2); z = coor * bz';

  [iswarp, evec] = mfs_warpck(x, y, z, geom);

# Warping correction

  if (iswarp)
     uv = [1, 2, 7, 8, 13, 14, 19, 20];
     r  = uv + 3;
     S  = speye(24, 24);
     Q  = mfs_warpQ(x, y, z);
     [R, Mz] = mfs_warpR(x, y, z, evec);
     R = repmat(R, 4, 1);
     S(3 : 6 : 21, uv) = [Q; -Q];
     S([1,  7,  8, 14], r) =  R;
     S([2, 13, 19, 20], r) = -R;
     S( 6,  4 :  5) = Mz(1 : 2);
     S(12, 10 : 11) = Mz(3 : 4);
     S(18, 16 : 17) = Mz(5 : 6);
     S(24, 22 : 23) = Mz(7 : 8);
     TE = S' * TE;
  end

end
