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

# usage: [coore, TE, T33] = mfs_esys(coor, opt)
#
# Input  coor(4, 3)   Coordinates of element nodes
#        opt          if 1, return 3 x 24 transformation matrix
#                     if 2, return 4 x 24 transformation matrix
#        geom         Geometrical data
#
# Output coore(4, 2)  Coordinates of element nodes in 2d element
#                     coordinate system
#        TE(:, 24)    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

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

# Check for warping

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

# Transformation matrix

  if (opt == 2 || iswarp)
     ix   = 1 : 3; incx = 4;  % Translations
     ir   = 4;                % Rotation
     TE1  = T33;
     TE   = sparse(16, 24);
  else
     ix   = 1 : 2; incx = 3;  % Translations
     ir   = 3;                % Rotation
     TE1  = T33(1 : 2, :);
     TE   = sparse(12, 24);
  end 

  iy = 1 : 3;

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

# Warping correction

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

end
