## 1 Introduction  
The numerical discrete element method (DEM) explicitly represents materials as assemblies of discrete particles or elements whose interactions dictate the overall behavior of the system [1]. In the DEM, each particle is assigned mass, geometry (which usually is as simple as a sphere [2] or a more complex shape [3]), and material properties, which are related to the contact interaction between particles governed by several contact laws. At the simplest level, the interaction is often modeled using an elastic spring, accompanied by tangential (shear) interaction with friction law [4]. To model solid materials, particles adhere to one another due to cohesive forces, assigning a cohesive bond between particles at the contact points. This approach is particularly suited for materials whose mechanical response is governed by the movement, contact, and breakage of individual grains or aggregates, such as granular soils [5], rocks [6], concrete [7], composites [8], and many other materials. Cohesive interactions enable crack initiation and propagation simulation by defining bonds that resist tensile and shear forces until failure [9]. When the force across a bond exceeds a critical value, the bond is broken, thereby simulating the initiation of microdamage. The accumulation of such breakages can eventually lead to macroscopic fracture. Recently, such an approach of modeling damage evolution by using DEM has been utilized widely [10]. Concrete fracture phenomena and their statistically varying character were investigated in [11], examining the complex processes of failure mechanisms, crack propagation, and damage evolution using the 2D discrete element model. In [12], a three-dimensional DEM particle-based cohesive crack model was developed to model the mixed-mode fracture process of brittle materials, aiming to simulate the material transition from a solid phase to a particulate phase. Damage and fatigue crack growth simulation of quasi-brittle materials were performed using a continuum damage model implemented in a discrete element code [13]. In [14], shale disk samples were studied under Brazilian test conditions using a 3D digital image correlation system (DIC) together with an inherently anisotropic model established within a DEM framework. Damage evolution and crack propagation of more advanced material, a fiber composite, was simulated in [8] concerning various fracture modes, such as fiber-matrix debonding, ply-to-ply delamination, matrix cracking, and fiber rupture. Thus, various DEM models are not limited to the analysis of granular materials, but can be effectively applied to studying cohesive powders [15], sea ice failure [16], the hardening/softening behavior of soil [17], starch agglomerates [18], high temperature resistant concrete [19], and coupled multi-physical problems [20]. Despite the efficiency of the mentioned DEM models in predicting the damage evolution of brittle materials, they are formulated based on regular/spherical discrete elements, which do not have any physical representations and were utilized only as the meshing structural feature. In the last several years, we can find more sophisticated approaches aimed at representing the damage of materials consisting of microstructural objects, such as grains [21,22] or aggregates [23,24,25,26]. The micro objects can be modeled as Voronoi polyhedra designed to resemble grains [21], accounting for the grain-scale heterogeneity observed in natural granular rocks. Grain-scale heterogeneity reproduced by Voronoi tessellation has also been discussed in [22]. Most of these grain-scale models exploit several reconstruction techniques to obtain the actual representation of the material microstructure, starting from optical microscopy [21,27] and ending with X-ray micro-computed tomography (micro-CT) imaging [23,24,26]. The exact arrangement of structural features, such as grains, particles, fibers, pores, or additional phases is crucial in the context of microcracking initialization and subsequent propagation.