This paper explores a modelling approach for two-dimensional (2D) and three-dimensional (3D) heterogeneous rock sample containing different types of inclusions base on the level set method (LSM). Taking the elliptical inclusion as an example, its distribution is considered as periodical and random, and the dimension is the same or random in the 2D finite element mesh. The interface between the matrix and inclusions is described by the functions suggested in this paper. Meanwhile, the modelling method for irregular inclusions is also provided using LSM for both 2D and 3D finite element mesh. In the 3D modelling, inclusions are assumed to be spheroid and the distribution also can be periodical or random. The main advantage of this modelling method is that the finite element mesh can remain the same without considering material distributions. In other words, the features of material distribution are separated from the mesh. But it should be pointed out that this method requires high computational resource. Finally, using LSM, the smeared crack model can be used to simulate the features of hydraulic fracture propagation in the heterogeneous rock materials with the same mesh.