We present a nondestructive approach to map the heterogeneous viscoelastic moduli from time harmonic motion via a constrained optimization strategy under the framework of finite element techniques. The adjoint equations are carefully derived to determine the gradient of the objective function with respect to the viscoelastic moduli. The feasibility of this inverse scheme is tested with simulated experiments under various driving frequencies. We observe that the overall strategy results in well-reconstructed moduli. For low frequencies, however, the mapped loss modulus is of inferior quality. To explain this observation, we analyze two simple one-dimensional (1D) models theoretically. The analysis reveals that the known displacement amplitude is less sensitive to the loss modulus value at low frequencies. Thus, we conclude that the inverse method is incapable of finding a well-reconstructed loss modulus distribution for low driving frequencies in the presence of noisy data. Overall, the inverse algorithms presented in this work are highly robust to map the storage and loss modulus with high accuracy given that a proper range of frequencies are utilized.