How do I rotate a face around an axis in normal orientation in Python?
After some reading I managed a solution, but I'm not sure how robust it is:
tangent = face.calc_tangent_edge_pair()
rot_axis = -face.normal.cross(tangent)

The tangent is one of the other normal axes and from that the last axis is calculated as its cross product with the face's normal.
This seems dependent on which edges are the longest (the documentation for BMFace.calc_tangent_edge_pair says: &quot;Return face tangent based on the two longest disconnected edges.&quot;), so I imagine there's a better way to approach this!