I am using the vector cross product to calculate each surface area and thus the outward normal. However, this approach fails when the surfaces are non-planar. So what is the best approach to calculate surface area and normal for a non-planar face of a hexahedron?
If the surface is nonplanar then the surface normal is a (hopefully continuous and differentiable) function of location. You will have to integrate the surface function S(x,y,z) over any two of the dimensions. If this is an engineering problem, then of course use a finite-grid approach and numerical integration. Then treat each grid element as a planar surface.