## 2 Overview of the LDEM Approach  
The LDEM approach modeled solids as a set of mass-less uni-axial elements and nodal masses, forming cubic cells of 20 bars and 9 nodes, as depicted in Figure 1a. Each node had three degrees-of-freedom, corresponding to its displacements in each coordinate direction. The lengths of each longitudinal and diagonal element were denoted as \(L_{n}\) and \(L_{d}=\sqrt{3}L_{n}/2\), respectively. The equations relating the properties of the LDEM elements to the elastic isotropic medium were as follows:  
\[\delta=\frac{9\
u}{4-8\
u},\qquad\alpha=\frac{9+9\delta}{2(9+12\delta)}, \tag{1}\]  
\[A_{n}=\alpha{L_{n}}^{2},\qquad A_{d}=\frac{2}{\sqrt{3}}\delta{\alpha{L_{n}}^{ 2}}, \tag{2}\]where \(E\) and \(\
u\) denote the [PERSON]'s modulus and [PERSON]'s ratio, respectively, of the material, whereas the initial axial stiffness values of the normal and diagonal elements are \(EA_{n}/L_{n}\) and \(EA_{d}/L_{d}\). Equations (1) and (2) were proposed by [PERSON] and [PERSON] in [26], and their derivations within the context of the lattice element method used here were complemented in [42]. For \(\
u=0.25\), the cubic array in Figure 1 yielded the exact representation of the isotropic continuum, while the other values of \(\
u\) caused small differences to appear in the shear values [26]. The motion equation was enforced at every node to obtain the following,  
\[\mathbf{M}\ddot{\mathbf{u}}+\mathbf{C}\dot{\mathbf{u}}+\mathbf{f}-\mathbf{p}=\mathbf{0}, \tag{3}\]  
where \(\mathbf{u}\) is the vector of the generalized nodal displacements, \(\mathbf{M}\) is the diagonal mass matrix, \(\mathbf{C}\) is the damping matrix, and \(\mathbf{f}\) and \(\mathbf{p}\) represent the internal and external forces acting on the nodal masses, respectively. If \(\mathbf{M}\) and \(\mathbf{C}\) are diagonal, the system is not coupled, and the explicit central finite differences could then be used to integrate the equations of motion in the time domain. As a result, since the nodal coordinates were updated at every time-step, large displacements were innately handled. The softening law for quasi-brittle materials proposed by [43] was implemented by the triangular-element constitutive relationship (OAB) of the LDEM longitudinal and diagonal elements (Figure 1b), which accounted for the irreversible effects of the crack nucleation and propagation. To ensure that the overall formulation was diagonal, the model's damping was assumed proportional to the mass, as follows:  
\[\mathbf{C}=D_{f}\mathbf{M}, \tag{4}\]  
where \(D_{f}\) is a damping coefficient for filtering unwanted higher frequencies in the numerical simulations. The most effective approach for its determination remained an open issue in LDEM applications. For example, in [44; 45], this coefficient was calibrated by the direct comparison with experimental measurements from global damping specimens. The numeric model was simulated under the same loading and boundary conditions from the experiments, and \(D_{f}\) was chosen so the model's logarithmic decrement matched the experimental values. An alternative approach was to modify the LDEM's uni-axial constitutive law by including a visco-elastic term that depended on the strain velocity. Here, \(D_{f}\approx 2\pi f_{m}\xi_{m}\) was estimated from the model's energy spectrum, where \(f_{m}\) and \(\xi_{m}\) represent the natural frequency and damping ratio of the highest energy peak, respectively. This approach was also used in [46], which explores the influence of \(D_{f}\) in more depth in the application of LDEM for modeling wave propagation in soils. Figure 1b describes the constitutive law assumed by the LDEM model for the force exerted by a bar as a function of its strain, where \(\varepsilon_{p}\) denotes the elastic limit strain, \(\varepsilon_{m}\) is the maximum strain experienced by the element up to that point, and \(\varepsilon_{r}\) is the critical failure strain.