This paper presents a method for reconstruction of maximum shear stress and stress trajectories from discrete data on principal directions. The domain is divided into smaller subdomains where stress potentials are assumed to be linear holomorphic functions. The functions obey continuity along element interfaces, which is used to form the first group of equations. The known data on principal directions are used in the second group of equations. Therefore, no stress magnitudes are involved in formulations, which eventually leads to a homogeneous system of linear algebraic equations. In order to make the system inhomogeneous an extra equation is added. It represents mean value of maximum shear stress over the domain. The reconstructed maximum shear stress, therefore, includes an arbitrary positive multiplicative parameter.