function [x, y, TE, T33] = mfs_esys(coor, opt, geom)

# usage: [x, y, TE, T33] = mfs_esys(coor, opt, geom)
#
# Input  coor(4, 3)   Coordinates of element nodes
#        opt          if 1, return 2 x 12 transformation matrix
#                     if 2, return 3 x 12 transformation matrix
#        geom         Geometrical data
#
# Output x(4)         x-coordinates of element nodes
#        y(4)         y-coordinates of element nodes
#        TE(:, 12)    Matrix to transform displacements 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 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

  x = coor * bx'; y = coor * by'; z = coor * bz';

# Check for warping

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

# Transformation matrix

  if (opt == 2 || iswarp)
     ix   = 1 : 3; incx = 3;
     TE1  = T33;
     TE   = sparse(12, 12);
  else
     ix   = 1 : 2; incx = 2;
     TE1  = T33(1 : 2, :);
     TE   = sparse(8, 12);
  end 

  iy = 1 : 3;

  for n = 1 : 4
      TE(ix, iy) = TE1;
      ix += incx; iy += 3;
  end

# Warping correction

  if (iswarp)
     Q  = mfs_warpQ(x, y, z);
     uv = [1, 2, 4, 5, 7, 8, 10, 11];
     if (opt == 1)
        S = sparse(12, 8);
        S(uv, :) = speye(8);
        S(3 : 3 : 12, :) = [Q; -Q];
     else
        S = speye(12);
        S(3 : 3 : 12, uv) = [Q; -Q];
     end
     TE = S' * TE;
  end

end
