Rock failure is intrinsically intertwined with heterogeneity and randomness. Although explicit and implicit numerical methods (e.g., image-based, Voronoi, and statistical heterogeneous models) have been developed to consider rock heterogeneity, it is still challenging to characterise rock heterogeneity and model rock heterogeneous fracture. This paper developed an alternative method to model rock heterogeneity and fracture. A grid nanoindentation test is carried out to investigate the heterogeneous distribution of micromechanical properties of granite. Then, a random field model is developed based on the nanoindentation test results. The mean value, coefficient of variance and correlation length of the random field for the granite specimen are determined. Further, the random field model is embedded into the combined Finite and Discrete Element Method (FDEM). Heterogeneous numerical models for Semi-Circular Bend (SCB) tests with pre-existing cracks are developed to simulate mode-I and mixed-mode fractures of rock. SCB experiments are also carried out to verify the developed method. Finally, the effects of random field parameters on rock fracture are investigated. The results show that Young’s modulus of granite from nanoindentation tests is highly heterogeneous and spatial correlated, and the heterogeneity can be characterised by a Gaussian random field with mean value 65.64 GPa, coefficient of variance 0.18 and correlation length 0.57 mm. The peak loads and crack trajectories from the developed numerical results agree well with those from the experimental results. The propagating crack is prone to shielding strong zones and passing weak zones in rock. The larger the correlation length and coefficient of variation are, the larger the variation of peak load is, indicating that the larger the size and the dispersion level of minerals tend to cause the higher uncertainty of the rock strength.