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

# usage: mfs_warpck(x, y, z)
#
# Input  x(4), y(4), z(4)   Nodal point coordinates in element
#                           coordinate system
#        geom               Geometrical data
#
# Output iswarp             0 if element is not warped
#                           1 if element is warped
#        evec(3, 4)         Vectors along element edges (optional)
#
# The function checks if a quadrilateral element is warped, and
# returns the vectors along the element edges.
#
# -------------------------------------------------------------------

  ib = [2, 3, 4, 1];  % Indices of edge end nodes

  if (isfield(geom, "wtol"))
     wtol = geom.wtol;
  else
     wtol = 1e-3;
  end

# Edge vectors

  coor = [x, y, z]';
  evec = coor(:, ib) - coor;
  ln   = norm(evec, "cols");

# Check warping

  lmin   = min(ln);
  iswarp = abs(z(1)) > wtol * lmin;

end
