•A detailed FDEM numerical method to simulate mechanical and fracturing responses of heterogeneous geomaterials with irregular inclusions is systematically developed.•A computational geometry method named CWSVM is proposed to control mesh quantity and quality.•A signed-distance-field-based discrete element method (SDF-DEM) is employed to approach the natural allocation and orientation of inclusions.•A combined constitutive model is proposed to consider the shearing hardening behaviour for the cohesive elements.•Effects of the interface strength on the mechanical and fracturing behaviours of inclusion-containing geomaterials are extensively discussed. In this paper, a detailed FDEM approach to simulate the mechanical and fracturing responses of heterogeneous geomaterials with irregular inclusions is systematically developed. The inclusion surface morphology is first obtained through 3D scanning techniques. A computational geometry method, the curvature-weighted sphere Voronoi method (CWSVM), is adopted to control the mesh quantity and quality and ensure the efficiency and accuracy of the FDEM numerical model. A signed-distance-field-based discrete element method (SDF-DEM) is employed to approximate the natural distribution and orientation of inclusions. Heterogeneous geomaterials with large inclusion contents (such as 60% and 70%) are generated effectively and efficiently through this approach. Next, to model the fracturing process, a finite discrete element method (FDEM) model is developed by integrating cohesive elements into the mesh in a fast and efficient manner. In addition, a combined constitutive model is proposed to consider the shear-hardening behaviour of the cohesive elements. The proposed numerical approach is verified through comparison with experimental results, including the shape of inclusions and mechanical responses of geomaterials. The results demonstrate that both satisfactory precision and low calculation costs can be achieved using the proposed algorithm. The consequent simulation performance is verified through comparisons of observations and numerical results with experimental results for failure patterns and mechanical behaviours. In addition, the effects of the strength of the interfaces between the inclusions and matrix on the mechanical and fracturing characteristics of inclusion-containing geomaterials are analysed quantitatively. The mechanical strength decreases rather than increases with increasing content of inclusions when the interface strength is less than the matrix strength.