## 1 Introduction  
Great engineering practices include the prediction and prevention of the propagation and initiation of cracks in structures with complex initial imperfections, such as rock masses with joints and concrete structures with early age cracks. When these structures are subjected to complex loading conditions, the existing cracks may propagate further, and new cracks may emerge. For these structures, analytical and empirical analyses are not sufficient to ensure their safety and durability. Numerical tools are highly preferable. With the understanding of continuum-discontinuous theory and developments in computing power, many sophisticated numerical methods have been proposed in recent decades. These methods can be built in a continuum or discrete framework [1; 2; 3] and can introduce damage degrees [4; 5] or crack openings (crack widths) [6; 7] as the new degrees of freedom. They can localize the damage [8; 9], consider nonlocal effects [10; 11], or even assume long-range forces [12; 13]. A crack can be explicitly represented by moving boundaries [14; 15; 16; 17] or implicitly embedded to avoid remeshing [18; 19; 20; 21]. The cracked domain can be discretized with elements [6; 22; 23] or particles [24; 25; 26]. In summary, these methods have advantages as well as disadvantages for different problems. Hence, a problem-oriented selection procedure would be helpful. Returning to our problem, to analyze the damaging processes of structures with initial imperfections, we need a method that is capable of capturing the initiation as well as the propagation of cracks. Since we do not focus on the stress state, such as the stress intensity factor around specific crack tips, numerical methods using remeshing [15; 27; 28], and nodal enrichment methods, such as the extended finite element method (XFEM) [18; 29; 30; 31] and the numerical manifold method (NMM) [32; 33; 34; 35], are not considered for simplicity. Furthermore, crack openings are very important for analyzing the durability of structures, and quasi-static loading conditions are considered. Hence, we do not use damage degree-based methods such as the phase field method [36; 37; 38; 39], mixed mode model [40; 41; 42; 43], equivalent lattice models [44; 45; 46; 47], and peridynamic-based methods [48; 49; 50; 51]. Finally, we hope that multiple cracks can be efficiently tracked simultaneously and that complicated crack tracking strategies [52; 53] can be avoided. The cracking elements method (CEM) [54; 55; 56; 57] is the chosen numerical tool. The CEM is a novel numerical approach belonging to the family of strong discontinuity embedded approaches (SDAs) [58; 59; 60; 61]. CEM uses Satically Optimal Symmetric (SOS) formulation [62]. Unlike the other SDAs, it introduces disconnected piecewise cracking segments appearing in the center point of each cracked element to represent crack paths, similar to the cracking particle method (CPM) [63; 64; 65; 66]. Hence, it does not need to distinguish crack tip elements and crack passing elements, which naturally captures the initiation and propagation of cracks. The crack orientation is determined locally, greatly reducing the computing effort. Moreover, [56] shows that a cracking element can be treated as a special type of finite element, which is formally like a 9-node quadrilateral element (or 7-node triangular element). Hence, this idea can be implemented easily in the conventional FEM framework. In this work, the CEM is used for simulating the damaging processes of brittle structures with joints and fissures. The numerically obtained results are compared with the experimental results, where the cracks do not always propagate from the tips of joints. In particular, uniaxial compression tests are considered, while we only use tensile damage criteria in our numerical model. On the one hand, satisfactory results are obtained in most cases where the cracking patterns observed in the experiments are numerically-reproduced.