I am working on an FEA mesh of hexahedron elements. The elemental level calculations involve finding the surface normals and area for each surface of a hex element. I preferred the vector cross product approach. But that doesn't seem to work when the faces are non-planar. Will the isoparametric transformation of each surface be a better option? Since I am already in a finite element setting.
If your vertices aren't coplanar, then your bounding surfaces can't be planes, and your volume is actually "worse" than an irregular hexahedron. Indeed, it isn't actually a polyhedron at all. At this point, having four vertices to a face stops being an over-specification of the surface normals, and starts being too little, in that there are an infinite number of curved surfaces passing through the four points.
This is where you need to take a step backwards and look at what caused you to ask the question in the first place. Since you're doing finite element analysis, it's very likely that the isoparametric trilinear basis through each vertex is actually the definition of the shape you're working with. As such, it's probable that the normals and surface areas you want are the ones arising for that object.
As a bonus, that probably means you can calculate the surface normals (in physical space) at locations given by parametric coordinates (the quadrature points for your method) rather than by physical coordinates. That would save you an inversion for every point on the surface where you need the normal (note again that since the surface isn't planar, its normal isn't constant).