Automatic hexahedral mesh generation plays a fundamental role in CAD/CAE fields, especially for IGA (isogeometric analysis). Recently, an automatic hex-mesh generation method has been proposed, which is based on the equivalence among three key concepts: colorable quadrilateral meshes, finite measured foliations and Strebel differentials.This work focuses on the computational aspect of Strebel differentials on high genus surfaces, which is based on graph-valued harmonic mapping. The algorithmic pipeline is as follows: first, the user inputs an admissible curve system, which induces a cylindric-decomposition graph; then the user specifies the lengths of the edges of the graph; third, the algorithm finds the unique harmonic map from the surface to the metric graph by a non-linear heat flow; finally, the harmonic map induces a Strebel differential, which can be further utilized to generate the hex-mesh.The method has solid theoretic foundation, which guarantees the existence and the uniqueness of the graph-valued harmonic map. The algorithm is capable of handling surfaces with complicated topologies, and producing all possible Strebel differentials on the surface. The user has full control of the combinatorial type and the geometry of the Strebel differential. The computational pipeline is automatic.The experimental results show the efficiency and efficacy of the algorithm, and demonstrate the great potential for automating the structured hexahedral mesh generation.