@pimms.value<EOL><INDENT>def face_normals(face_coordinates):<DEDENT>
u01 = face_coordinates[<NUM_LIT:1>] - face_coordinates[<NUM_LIT:0>]<EOL>u02 = face_coordinates[<NUM_LIT:2>] - face_coordinates[<NUM_LIT:0>]<EOL>if len(u01) == <NUM_LIT:2>:<EOL><INDENT>zz = np.zeros((<NUM_LIT:1>,tmp.shape[<NUM_LIT:1>]))<EOL>u01 = np.concatenate((u01,zz))<EOL>u02 = np.concatenate((u02,zz))<EOL><DEDENT>xp = np.cross(u01, u02, axisa=<NUM_LIT:0>, axisb=<NUM_LIT:0>).T<EOL>norms = np.sqrt(np.sum(xp**<NUM_LIT:2>, axis=<NUM_LIT:0>))<EOL>wz = np.isclose(norms, <NUM_LIT:0>)<EOL>return pimms.imm_array(xp * (np.logical_not(wz) / (norms + wz)))<EOL>
mesh.face_normals is the (3 x m) array of the outward-facing normal vectors of each
  triangle in the given mesh. If mesh is a 2D mesh, these are all either [0,0,1] or
  [0,0,-1].