Reaction diffusion is the most common growth modelling methodology due to its simplicity and consistency with the biological tumor growth process. However, current extensions of the reaction diffusion model lack one or more of the following: efficient inclusion of treatments' effects, taking into account the viscoelasticity of brain tissues, and guaranteed stability of the numerical solution. We propose a new model to overcome the aforementioned drawbacks. Guided by directional information derived from diffusion tensor imaging, our model relates tissue heterogeneity with the absorption of the chemotherapy, adopts the linear-quadratic term to simulate the radiotherapy effect, employs Maxwell-Weichert model to incorporate brain viscoelasticity, and ensures the stability of the numerical solution. The performance is verified through experiments on synthetic and real MR images. Experiments on 9 MR datasets of patients with low grade gliomas undergoing surgery with different treatment regimens are carried out and validated using Jaccard score and Dice coefficient. The growth simulation accuracies of the proposed model are in ranges of [0.673 0.822] and [0.805 0.902] for Jaccard scores and Dice coefficients, respectively. The accuracies decrease up to 4% and 2.4% when ignoring treatment effects and the tensor information, while brain viscoelasticity has no significant impact on the accuracies.