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

# usage: [coore, TE, T33] = mfs_esys(coor, geom)
#
# Input  coor(3, 3)   Coordinates of element nodes
#        geom         Geometrical data
#
# Output coore(3, 2)  Coordinates of element nodes in 2d element
#                     coordinate system
#        TE(18, 18)   Matrix to transform displacements and rotations
#                     to element coordinate system
#                     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 triangular shell elements.
#
# --------------------------------------------------------------------

# Basis vectors of element coordinate system

  bx = coor(2, :) - coor(1, :);
  bx = bx / norm(bx);
  c2 = coor(3, :) - coor(1, :);
  bz = cross(bx, c2);
  bz = bz / norm(bz);
  by = cross(bz, bx);

  T33 = [bx; by; bz];

# Coordinates in element coordinate system

  relcoor = coor - coor(1, :);
  coore   = relcoor * [bx', by'];

# Transformation matrix

  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 : 3
      TE(ix, ix) = T66;
      ix += 6;
  end

end
