{'fraction_non_alphanumeric': 0.03327660609007476, 'fraction_numerical': 0.009967787029721023, 'mean_word_length': 4.7109510586601875, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 9, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Comparative study of multiscale computational strategies for materials with discrete microstructures

The evolution of local defects such as dislocations and cracks often determines the performance of engineering materials. For a proper description and understanding of these phenomena, one needs to descend to a very small scale, at which the discreteness of the material emerges. Fully-resolved discrete numerical models often suffer from excessive computing expenses when used for application-scale considerations. More efficient multiscale simulation procedures are thus called for, capable of capturing the most significant microscopic phenomena while being computationally tractable. Two broad classes of methods are available in the literature. The first class considers the fully-resolved discrete system, which is subsequently reduced through suitable mathematical tools such as projection and reduced integration. The second class first homogenizes the discrete system into an equivalent continuum formulation, into which the main phenomena are added through specific enrichments. This paper provides a thorough comparison of the two different modeling philosophies in terms of their theory, accuracy, and performance. To this goal, two typical representatives are adopted: the Quasicontinuum method for the first class, and an effective continuum with an embedded cohesive zone model for the second class. Two examples are employed to demonstrate capabilities and limitations of both approaches: dislocation propagation and pile-up against a coherent phase boundary, and a three-point bending test of a concrete specimen with crack propagation. In both cases, the accuracy of the two methods is compared against the fully-resolved discrete reference model. It is shown that whereas continuum models with embedded cohesive zones offer good performance to accuracy ratios, they might fail to capture unexpected mechanical behavior. The Quasicontinuum method offers more flexibility at a higher computational cost.


Introduction
Numerical modeling of material behavior is an important ingredient in engineering. At sufficiently large scales, most materials can be considered as a continuum, whereas at smaller scales, materials such as foams, textiles, concrete, or paper, reveal discrete microstructures. At the smallest scale, i.e., nano-scale, all materials become discrete and their behavior is dictated by the underlying atomic lattice. When modeling of localized mechanical phenomena such as for instance the nucleation and propagation of dislocations, cracks or other defects, is of interest, the inherently non-local behavior of the underlying discrete microstructure comes into play. To capture those non-localities, accurate numerical models may be required.
The most accurate option is to use discrete models that fully resolve the microstructure and hence inherently incorporate the underlying non-locality at the level of the individual interactions. Typical examples at different scales are molecular statics [65,20,55] or dynamics [26,68,79] for atomic systems, or lattice or beam networks for mesoscopically discrete materials. Sophisticated discrete material models for various materials such as paper [43,19], textile [13], fiber-reinforced composites [6,80], other fibrous materials [77,59], or concrete [63,25,42] can be found in the literature. For typical engineering applications, however, full-scale models suffer from prohibitive computational costs requiring simplifications, especially in problems with a large separation of scales between the application length scale and the length scale of the underlying lattice, or in problems in which multiple model evaluations are required such as optimization, parameter identification, or stochastic modeling.
Two broad classes of conceptually distinct methodologies are available to provide computationally efficient numerical tools for such cases. The first class relies on considering directly the underlying model, which is subsequently reduced by suitable numerical and mathematical techniques such as projection on a reduced kinematical basis and reduced integration. The second class of methods first derives equivalent homogeneous governing equations from the underlying discrete microstructure. Afterwards, suitable enrichments are added to model the key mechanical phenomena. The resulting system is then solved using standard numerical techniques for continuum problems such as the Finite Element (FE) method. For the purposes of this paper, the QuasiContinuum (QC) method, originally introduced by Tadmor et al. [71], is adopted as a typical representative for the first class of approaches, whereas a continuum formulation with an embedded cohesive zone model is considered as a representative for the second class of approaches. The reduction of both classes is implemented at different levels of the problem formulation, as depicted in Fig. 1. In the first approach, the original system is reduced through suitable numerical and mathematical techniques, solving subsequently a simplified system (the left-bottom path). In the second approach, the original system is homogenized first to provide an equivalent continuum model with embedded key features, which is subsequently solved using standard techniques for continuum problems (the top-right path). (b) A graphical representation of the individual steps required by the two simplification approaches.

Quasicontinuum based approach
In the QC approximation, the original complexity of the full underlying problem is reduced to simplified model in two steps: (i) Interpolation is used to define the positions of all atoms or lattice sites of the original system from a set of representative atoms, the so-called repatoms; (ii) Summation efficiently estimates the stored energy and internal forces of the system in which the interpolation step has been applied, based on a set of the so-called sampling atoms.
The interpolation step considers only a small set of repatoms to represent the kinematic behavior of the entire system, typically chosen as vertices of an overlaid triangular Finite Element (FE) mesh. To accurately capture any localized phenomena, the considered domain is split into two parts. The first part is the area of high interest where all atoms of the underlying lattice are considered as repatoms and thus also as nodes of the triangular mesh, i.e., the overlaid triangulation is fully refined down to the level of the underlying lattice. Elsewhere, interpolation is adopted to coarse-grain the underlying lattice; in this region individual triangular elements contain multiple atoms or lattice sites. The fully-resolved region may be adaptively adjusted during the simulation to capture the relevant evolution of the microstructure. The main advantage of the QC approach is that since the response of the reduced model is governed by the underlying 'exact' discrete model, one expects a high fidelity, which is limited only by the coarsening away from the region of interest. The price to pay is (possibly) a more expensive simulation.
The QC method was originally introduced in 1996 [71,72] for atomic systems with longrange conservative interaction potentials and has been successfully used for the simulation of localized phenomena, such as crack nucleation and propagation [37,52], motion and interaction of dislocations [57], or processes related to nanoindentation [38,69]. Later, the QC method was extended to regular discrete lattice networks with short-range nearestneighbor interactions with both conservative [8] and non-conservative [10,11] interaction potentials including dissipation and fibre sliding. A further extension was provided for planar beam lattices [12,7] and periodic beam lattices in a co-rotational framework [58]. Additional research focused on goal-oriented adaptivity [2,46], a meshless QC method [40], or an energy-based variational formulation for regular lattices with plasticity [60] and localized damage [62,61]. Recently, the QC method was also extended to irregular lattices [50] or polymer networks [31]. A generalization to metallic lattice materials was introduced in [21], where different finite element shape functions are used for different types of lattice nodes.

Homogenization based approach
In the second approach, homogenization of the underlying discrete system is first carried out, from which continuum governing equations follow. Usually a Cauchy continuum is adopted, which neglects any non-local interactions. A specific enrichment thus needs to be added, which is usually implemented through an appropriate interface model, e.g., a cohesive zone model for cracks. These interfaces allow for mutual displacements in tangential and normal directions, governed by appropriate effective constitutive laws.
The positioning and orientation of such interfaces is assumed to be known a priori and fixed hereafter, although more general models can be found in the literature, involving cohesive zones in between all elements [78], or an extended FE method with cohesive zone formulation [14,4,76]. Often, the interfaces positions can be set a priori on the basis of available knowledge of the underlying physics and to limit the computational expense. A crack path can in certain cases be estimated with sufficiently high accuracy to be fixed in all simulations, for instance. Thus, by design, this approach is capable of capturing only those phenomena that are considered in the homogenized model. In situations in which the positions of interfaces cannot be estimated in advance, homogenization based approaches are usually no longer advantageous. Even though generalized homogenized models (which, e.g., consider cohesive zones in between all elements) can be used even without any prior information on the position of the expected localized phenomena, it is usually computationally less expensive directly employ adaptive QC method in such situations, since induced computing costs might be comparable or exceed those of the QC method or even of a full solution. Although many physical phenomena are excluded from the homogenized modeling, one generally expects this approach to be more efficient. The material parameters of the homogenized system as well as properties of the cohesive interfaces are fully determined by the underlying lattice physics. In most cases, a linear elastic isotropic material model is sufficient to provide an accurate description of the lattice behavior far from any localized phenomena, whereas in the close vicinity of a crack or a dislocation large deformations and strains may occur; these are, however, neglected in this contribution for simplicity. Interface constitutive laws range from simple closed-form uncoupled to numerical and coupled laws derived directly from numerical homogenization. Typical representatives of the homogenization based approaches, both considered in this paper, are the Peierls-Nabarro model introduced in [56,54,64,34,15] and cohesive zone model discussed in [5,29,18,30,35,67]. Yet another representative is the Discrete Dislocation Dynamics (DDD) approach, e.g., [32,23], which is not considered here further.

Objectives
The goal of this paper is a thorough comparison of the above two modeling philosophies in terms of their theoretical framework, accuracy, and numerical efficiency. To this end, two examples are considered with microstructures at the nano-and meso-scale level. At the nanoscale, dislocation propagation and pile-up against a coherent phase boundary in a twodimensional atomic system with a hexagonal stacking of atoms is adopted. At the mesoscale, a three-point bending test of a concrete specimen with crack propagation through a twodimensional regular X-braced lattice structure is examined. In both cases, a Quasicontinuum model and a FE approach equipped with cohesive zone interfaces, using either the Peierls-Nabarro model or a traction-separation law, are compared against a fully-resolved discrete model. The accuracy and efficiency of both homogenization techniques are discussed and evaluated.

Outline and notation
The remainder of this paper is organized as follows. Section 2 considers the molecular statics problem of dislocation propagation, starting with the theory of the fully-resolved discrete atomic problem, followed by a Quasicontinuum formulation, and closing with the FE Peierls-Nabarro homogenization based approach. In Section 3, an analogous comparison is made for the case of a three-point bending test considered at the mesoscale, mimicking quasibrittle damage in a concrete specimen. The paper closes with a summary and conclusions in Section 4.
Throughout the paper, scalars are denoted a, vectors a, position vectors in a twodimensional continuum x = x e x + x e y , and atom position vectors r = r x e x + r y e y . Secondorder tensors are denoted σ, matrices A, column matrices a, scalar products of two vectors a · b = a i b j , single contractions of two second-order tensors σ · ε = σ ik ε kj e i e j , double contractions of two second-order tensors σ : ε = σ ij ε ji , whereas a transpose is denoted as σ T , where σ T ij = σ ji . The gradient and divergence operators are denoted as ∇ a( x) = ∂a j ( x) ∂x i e i e j and ∇ · a( x) = ∂a i ( x) ∂x i , respectively. Throughout this contribution, Einstein's summation convention is adopted on repeated indices i and j, and e x and e y denote the basis vectors of a two-dimensional Cartesian coordinate frame.

Full atomic model
To analyze a dislocation pile up against a material interface in two dimensions, a shear test with two different material phases is considered. Dislocation dipoles are nucleated at the center of the specimen, and are pushed towards boundaries by externally applied shear Inside Ω, hexagonal lattice structures of equal spacings a 0 are considered in both phases, comprising a total of n ato atoms α, collected in an index set N = N A ∪ N B where N A and N B are sets of atoms associated with the individual phases. The crystal orientation is chosen such that one set of glide planes is oriented perpendicular to the two vertical coherent phase boundaries, whereas the other sets of glide planes are inclined by an angle of ±60 • . The interatomic pair potentials employed within both phases, denoted φ A and φ B , are of the Lennard-Jones (LJ) type. To account for nearest and next-to-nearest interatomic interactions only, a cut-off radius r cut = 2.5 r m is adopted, as indicated in Figs. 2a and 3a by the dashed circles. To ensure smoothness of the potential at the cut-off radius, cf. Fig. 3b, the LJ potential is modified according to where r is a scalar measuring the distance between a pair of atoms, and is the standard LJ interatomic pair potential. In Eq. (2), ε denotes the depth of the energy well and r m the distance for which the interaction energy reaches its minimum, as indicated in Fig. 3b. For further details see, e.g., [70]. Mixing of the interatomic potentials across the phase boundary Γ pb is considered through averaging, i.e., where r α = r α x e x + r α y e y denotes a current position of an atom α, r αβ = r αβ 2 denotes the Euclidean distance between a pair of atoms α and β, r αβ = r β − r α is a vector of their relative positioning, and Material parameters are specified in Tab. 2, which are all normalized with respect to ε A and r A m . If only the nearest-neighbor interactions are considered and there are no surface effects within a homogeneous crystal, i.e., an infinite homogeneous crystal is assumed, the lattice (a) fully-resolved discrete atomic system    spacing a 0 = r m yields a stress-free configuration. For next-to-nearest interactions, however, such a lattice spacing invokes a pre-stress within the lattice, as calculated by a non-vanishing Virial stress tensor σ defined as where V denotes a volume considered in the deformed configuration over which the stress is computed, see, e.g., [70, Section 11.5.2] for more details. To attain a stress-free state, the system reduces its initial interatomic spacing a 0 = r m to a 0 = 0.99296702 r m . Although such a change is geometrically negligible, it has significant influence on the macroscopic mechanical behavior of the resulting system. For the application of the boundary conditions the behavior of an ideal, defect-free, and linear-elastic specimen is assumed under a state of constant shear stress τ , i.e., As illustrated in Fig. 2b, this is achieved by prescribing shear strains γ inversely proportional to the effective shear moduli µ associated with the individual phases. In particular, the shear angles correspond to tγ in Phase A, and tγ/ρ in Phase B, where ρ is a material contrast ratio, and γ is a target shear angle in Phase A, and t ∈ [0, T ] is a parametrization pseudotime. The requirement of constant shear stress translates to prescribing atom positions of the boundary atoms as where N BC ⊂ N are the atoms to which boundary conditions are applied, i.e., atoms positioned within the r cut distance from the specimen's boundary ∂Ω, indicated by the hollow circles in Fig. 2a.
The mechanical configuration of the system at each time step t k for a time discretization of the time horizon t ∈ [0, T ], 0 = t 0 < t 1 < · · · < t n T = T, is governed by the minimization of the total potential energy of the system, i.e., Here, r k = r(t k ) is a relaxed configuration at a time instant t k , r is a column matrix storing the components of the position vectors in the deformed configuration r α corresponding to all atoms α = 1, . . . , n ato , v denotes an admissible configuration of the system reflecting the boundary conditions of Eq. (7) applied at a time instant t k , and is the internal free energy 1 . No externally applied forces are considered, as only the prescribed displacements are used to load the atomic system, recall Eq. (7). The minimization problem of Eq. (9) is solved using a Trust-region algorithm, see, e.g., [24], which has been implemented within an in-house code.
To avoid potential occurrence of multiple glide planes, the numerical solver is initiated at the beginning of each time step t k towards the preferred glide plane Γ gp = { x ∈ R 2 : y = 0} by means of the analytical Volterra solution; see Fig. 2a, where the preferred glide plane is depicted by the dotted line. A more detailed description of the same initialization procedure can be found in [15, Section 2.2].

Quasicontinuum model
The Quasicontinuum (QC) methodology is a concurrent multiscale technique, originally introduced in [71]. The main idea consists in combining an accurate but expensive full atomic description inside regions of high interest with a cheaper continuum-like approximation elsewhere. This is achieved through a finite element mesh overlaid on top of the underlying lattice, as shown in Fig. 4, where all mesh nodes (vertices of linear triangular elements) coincide with atoms. This set of atoms, the so-called repatoms, determines positions of all the remaining atoms through interpolation In Eq. (11), r rep is a column storing positions of all repatoms, whereas Φ = [ϕ 1 , . . . , ϕ 2nrep ] stores the interpolation basis functions ϕ i associated with the triangular mesh, see [70] for more details. By gradually refining the triangulation, a seamless transition between the fully-resolved discrete description in the region of high interest and an efficient continuum description elsewhere is achieved. A common practice is to choose the mesh such that the number of repatoms n rep is much smaller than the total number of atoms n ato , reducing substantially the number of Degrees Of Freedom (DOFs). This in turn yields computational savings, since the minimization problem of Eq. (9) is carried out with respect to r rep instead of r. Repatoms are shown as black dots, sampling atoms as red dots, and all remaining atoms as gray dots. In the fullyresolved region the triangulation matches the underlying atomic structure (gray triangles), whereas elsewhere interpolation is used (blue triangles).
In order to construct the potential energy V in Eq. (10) all lattice sites need to be visited, which is expensive and unnecessary. The second Quasicontinuum step is thus introduced to approximate the energy by where the energy associated with a lattice site α is specified as and where w α is a summation weight associated with each sampling atom α ∈ N S . A set of all sampling atoms, N S (indicated in Fig. 4 as red dots), is selected to accurately represent the energy of the entire system. The summation rule is tailored to a particular interpolation scheme used, for which multiple options exist, as reported, e.g., in [1,12,44]. In what follows, the central summation rule of Beex et al. [9] is used. A QC method thus minimizes the approximate energy of Eq. (12) with respect to the reduced kinematic variable r rep related to the configuration of the entire atomic system through interpolation of Eq. (11). The area of high interest, i.e., the fully-resolved region and hence the associated triangulation, need to adaptively evolve in time to properly accommodate dislocation motion while retaining QC efficiency. To this end, various adaptive criteria for atomic lattices can be found in the literature, see, e.g., [53,58,46,73,41,3]. In this work, the Zienkiewicz-Zhu (ZZ) error estimator [81,82,66] is used as a refinement criterion, which is applied to the deformation gradient tensor computed within each element, see also [70,Section 12.6.5]. That is, the local error inside each element is estimated and compared against a pre-selected threshold value ZZ tr , to determine elements which need refinement. The fully-resolved region is updated sequentially as follows: (i) At each time increment, the QC system is equilibrated for a fixed fully-resolved region; (ii) The triangulation is checked by the mesh refinement criterion (using the ZZ error indicator) and refined if necessary; (iii) All required atoms are added as repatoms, the set of sampling atoms is amended, the interpolation matrix is updated, and the equilibrium is restored again.
This procedure is repeated until the mesh refinement criterion is satisfied for all mesh elements, proceeding subsequently to a new load increment in (i). Since in the here considered problem the same dislocation path is followed by trailing dislocations, no coarsening is considered for the atomic problem.

Peierls-Nabarro formulation
The homogenization based simplification of the full atomic problem of Section 2.1, shown in Fig. 2a, employs the FE Peierls-Nabarro model (abbreviated FE-PN; for more details see [17]). The problem domain Ω is split into two regions, Ω A and Ω B . These regions are separated by perfectly and fully coherent phase boundaries of zero thickness Γ pb , which are normal to e x . Displacement and traction continuity conditions are enforced on Γ pb as where σ is the Cauchy stress tensor, and u is the displacement field. Both phases Ω i , i ∈ {A, B}, are modeled as an elastic continuum enriched with a Peierls-Nabarro glide plane Γ gp . Any dislocation motion, and in fact all non-linearity in the system, is confined to this plane, inducing additional energy through a corresponding misfit potential. The glide plane is considered perpendicular to the phase boundary Γ pb , i.e., oriented along e x . Γ gp is assumed as continuous throughout the entire domain Ω, horizontally splitting each phase Ω i into two subdomains Ω i ± , i ∈ {A, B}, as shown in Fig. 5. A shear deformation is prescribed to the considered domain as a function of time t ∈ [0, T ] on the external boundary ∂Ω, inducing initially constant shear stress throughout the specimen according to Eq. (6).
In analogy to the full atomic system and the QC formulation, the evolution of the FE-PN model is governed by minimization of the total potential energy (specified per unit thickness because plane strain conditions are assumed), i.e., In Eq. (15), u k = u(t k ) is the relaxed displacement field at a time instant t k and v is a kinematically admissible displacement field respecting the kinematic constraints prescribed on ∂Ω. The total potential energy is expressed as where ψ e is the elastic strain energy density (considered inside Ω i ± ), ψ gp is the glide plane potential (localized along the glide plane Γ i gp ); any external forces have been neglected. Within all subdomains Ω i ± , small strain linear elasticity is assumed, i.e., where is the small strain tensor, and ∇ denotes the gradient operator in the reference configuration. The energy density of Eq. (17) yields the linear elastic constitutive law in which D is a phase-specific isotropic fourth-order elasticity tensor.

Calibration
The elasticity tensor D used in the FE-PN model in Eqs. (17) and (19), which is fully isotropic due to the considered hexagonal lattice, is obtained through homogenization of the atomic lattice as follows (see also [70]). First, a numerically homogenized stiffness tensor D atm , with components is computed on the basis of the periodic unit cell of the atomic lattice. In Eq. (20), V denotes the volume of the simulation cell in the deformed configuration, whereas r αβ i (i = 1, 2) are the components of the relative position vectors r αβ , cf. Eq. (5) and the discussion thereof.
For nearest and next-to-nearest interactions the stiffness tensor D atm has been derived analytically for a unit hexagonal lattice by assembling the contributions of all relevant atoms (within the dashed circle in Fig. 3a). The computed parameters (tensor components) for Phase A are listed in Tab. 3 (the corresponding atomic properties are specified in Tab. 2); the parameters of Phase B are obtained simply by scaling those of Phase A with the material contrast ratio ρ ∈ {1.4, 4.0}. Next, an elastic isotropic plane strain stiffness tensor D(E iso , ν iso ) is considered with the Young's modulus E iso and Poisson's ratio ν iso , in which effective parameters The glide plane potential ψ gp is a function of the disregistry profile ∆, i.e., of the displacement jump across the glide plane Γ gp , which splits into a tangential ∆ t and a normal ∆ n part, i.e., ∆ = ∆ t e t +∆ n e n . To capture the effect of periodicity of the underlying atomic lattice, ψ gp is a non-convex periodic function of ∆ t with a period a 0 . Different expressions have been introduced for ψ gp in the literature, as discussed and compared, e.g., in [15]. Here it is based on the Generalized Stacking Fault Energy (GSFE) of the underlying atomic lattice, see, e.g., [75] where the GSFE was obtained from atomic calculations, or [28] where density functional theory was employed. In this paper, the GSFE is constructed through atomic calculations using a rectangular simulation box of size 20a 0 ×12a 0 √ 3 with periodicity conditions between the vertical boundaries and free horizontal surfaces, as shown in Fig. 6a. A lattice with the stress-free spacing a 0 is considered, and the upper part is rigidly displaced as ∆ = ∆ t e t +∆ n e n . The GSFE computed for Phase A is shown in Fig. 6b, whereas the GSFE for Phase B is obtained again by multiplying with the corresponding material contrast ratio ρ. Within the specimen domain Ω, the glide plane potential is considered as a piecewise constant function. This is possible because the lattice considered in both phases is perfectly aligned, having the same lattice spacing a 0 . If the lattices within the different phases were not perfectly aligned with the same lattice spacing, or if they were oriented differently with respect to each other, additional simulations of the phase interface might be required.
Standard finite elements are used to discretize the elastic regions Ω i ± . The glide plane Γ gp is discretized by interface elements which are inserted between the bulk elasticity elements above and below it. Mechanical equilibrium is established by minimizing the total potential energy Ψ of Eq. (16) with respect to the DOFs of the considered FE triangulation. The non-convexity of ψ gp is addressed by a truncated Newton optimization algorithm elaborated in [16]. Individual dislocations are initialized in analogy to the approach used for the fullyresolved atomic system, as described at the end of Section 2.1; for more details see [17, Section 3.2].

Results and comparison
The behavior of the fully-resolved atomic system of Section 2.1 (referred to as Full) is described first, considering the material contrast ratio ρ = 1.4. Two QC systems of Section 2.2 are considered: (i) QC with a fixed mesh that is fully refined along the entire glide plane (referred to as QC fix), and (ii) QC with a mesh that has only a small fullyresolved region situated around the Frank-Read source equipped with adaptivity and a Zienkiewicz-Zhu error indicator of critical threshold ZZ tr = 0.001, referred to as QC ZZ. The FE-PN model of Section 2.3 is considered for only one discretization with a fixed mesh refined along the glide plane down to an element size h = a 0 /16. Although such an excessively small element size is not strictly necessary, see [17,15] where 4-8 times coarser mesh was employed to obtain adequate macroscopic results such as transmission stress, it is used here to test the best accuracy of the FE-PN method, in particular to accurately capture the shape of the dislocation core. For coarser meshes a substantially better speed-up can be expected compared to the results presented here. The initial triangulations associated with the individual systems are shown in Fig. 7, whereas the corresponding (initial/final) number of DOFs are listed in Tab. 4. All reduced methods are compared against the full atomic system in terms of the total potential energy evolution, dislocation positions, critical transmission stresses, disregistry profiles, and computational effort. Upon loading, for the fully-resolved atomic system at a certain time instance a critical resolved shear stress is reached (the activation stress of a Frank-Read source), resulting from the shear deformation applied on the remote boundary. At that moment, a new dislocation dipole is nucleated, which moves symmetrically under the increasing load and local Peach-Koehler force towards the phase boundary, where it is obstructed as a result of the phase contrast. By increasing the applied shear further, new dislocations are emitted and a dislocation pile-up is established in front of the interface. Further increasing the applied shear induces dislocation transmission into the neighboring Phase B. For visualization purposes, the quantity called Local Lattice Disregistry associated with an atom α (LLD α ), represents how much the current configuration of the hexagonal lattice around this atom differs from its initial configuration, and is defined as a sum of differences between the reference (undeformed) positions of its six nearest-neighbor atoms, r γ t , and their current deformed positions located closest to their initial positions, r β , i.e., where min ||r γ t − r β || 2 denotes the Euclidean distance of a one particular theoretical position from its closest atom β in the actual deformed configuration, and B α is the initial set of nearest-neighbors associated with an atom α. The LLD indicator of Eq. (23) is shown in color for all atoms of the full atomic system in Fig. 8, corresponding to the case of a dislocation pile up with subsequent dislocation transmission (i.e., ρ = 1.4, Figs. 8a and 8c). Global energy evolution profiles for all models are shown in Fig. 9a, where a good agreement is achieved, with the maximum energy error below 3% excluding jumps in energy evolutions. These jumps, present mostly in the early stages of the evolutions and observed mainly for the Full and both QC models, result from the initiation of new dislocations, which temporarily increases the energy. Upon further loading, the initially immobile but stable dislocations start to propagate towards the phase boundary and the energy drops. Similar behavior is not observed for the FE-PN model due to a vanishingly small Peierls-like barrier, as a consequence of a rather fine discretization (recall that h = a 0 /16 has been adopted on the glide plane).
The corresponding evolutions of dislocation positions along the glide plane Γ gp are shown in Fig. 9b against the externally applied normalized shear load τ /µ A . Here we notice that all reduced methods cluster around the results corresponding to the full solution, suggesting good agreement. If there is the same number of dislocations present in the compared models, then the maximum error in the dislocation position does not exceed 10 a 0 once the dislocations start propagating towards the interface. However, more significant discrepancy in the dislocation position can be observed at early stages after initiation. All dislocations in the initiation is very sensitive to stress distribution along the glide plane, differences in the Full and both QC methods are also non-negligible. Numerical settings, such as discretization and interpolation, might have an impact on the evolution of the initial dislocations, affecting their final trajectories. We further notice that the critical transmission stress is underestimated by more than 30% by the FE-PN method and slightly overestimated by the QC ZZ method. For clarity, the transmission stresses corresponding to the individual methods are listed in Tab. 4. The significant underestimation of the critical transmission stress by the FE-PN method originates from several reasons. The first reason is an inaccurate representation of the dislocation core structure by the FE-PN model (because of adopted linear elasticity neglecting large deformations, discreteness, and non-locality of the underlying lattice), which may play a dominant role for dislocation-interface interactions. Dislocation cores obtained for the Full and FE-PN models are compared in Fig. 10, where we clearly see that the dislocation core provided by the FE-PN model is almost symmetric with respect to the horizontal glide plane unlike the strongly asymmetric core obtained from the Full model. The second reason is a missing Peierls barrier causing more compressed pile-up, resulting in higher stresses acting on the leading dislocation in the FE-PN model. The last but probably the most important reason is that in the atomic models, the leading dislocation dissociates partially into the secondary inclined glide planes (associated with dislocation reflection, discussed later in this section, see also Fig. 8d). As a result, the energy of the dislocation is split between the primary and secondary glide planes. On the other hand, because in the FE-PN model the entire energy is constrained to the primary glide plane associated with the transmission, the dislocation is located closer to the interface and dislocation transmission thus happens earlier. Additional simulations (not presented here for brevity) confirm that the shear stress needed for dislocation transmission increases significantly once the secondary inclined glide planes are introduced into the FE-PN model.  Normalized tangential ∆ t /a 0 and normal ∆ n /a 0 disregistry profiles, expressed as a function of the normalized coordinate x/a 0 along the glide plane Γ gp , are shown in Fig. 11 for two load levels τ = 0.0068 µ A and τ = 0.01093 µ A . In all discrete models (Full, QC fix, and QC ZZ), the disregistry profiles are evaluated from the displacements of individual atoms along the glide plane Γ gp , whereas for the FE-CZ model the displacement discontinuity vector ∆ is plotted directly. From the presented results we conclude that, consistently with Fig. 9, at the lower applied shear level (τ = 0.0068 µ A ) the QC fix solution is very similar to the reference Full solution, whereas the FE-PN method achieves better accuracy as compared to the QC ZZ. A different situation occurs at later stages of loading (τ = 0.01093 µ A ), where the best accuracy is achieved by the QC ZZ method, while the FE-PN method leads to the highest error. This error is quantified in Tab. (4) by the Dislocation Position Error (DP E). For each reduced model, the DP E value represents how much on average the position of dislocation in this model differs (relative to a 0 ) from the exact solution. The averaging is realized for both distinct dislocations as well as distinct time steps, i.e., where x Full i d ,is and x • i d ,is denote horizontal position of i d -th dislocation at i s -th time step corresponding to the Full and one of the effective models (QC fix, QC ZZ, or FE-PN). n s denotes the number of time steps used for the error evaluation, and n d is the total number of dislocations at the given time step. In Tab. (4), the DP E quantity is evaluated as an average difference in position of three dislocations in the last two steps before transmission of the FE-PN model (τ = 0.01093 µ A ). The most accurate QC ZZ method exhibits an average difference of 5.83 a 0 , while the highest average difference of 8.92 a 0 can be observed for the FE-PN method.
The computational performances obtained for the individual models are summarized in Tab. 4 in terms of the number of DOFs and computing times. Both QC models achieve significant reductions in the number of DOFs, which also results in considerable computational savings. The QC fix requires less computing time compared to the QC ZZ method because in the latter a significant amount of time is spent on mesh refinement procedures. The FE-PN model, on the other hand, uses relatively fine discretization along the glide plane (recall that h max = a 0 /16 to capture dislocation core accurately), requiring even a higher number of DOFs as compared to the full atomic model. Only a five-fold speed-up is thus obtained, which can be significantly improved by employing much coarser meshes. Additional simula-tions show that 4-8 times coarser meshes relative to the fine h max = a 0 /16 discretization do not change the results by more than 5% while approximately a ten-fold speed-up is achieved, cf. [17,15].
Although the presented example was set to show dislocation pile up with subsequent dislocation transmission, another, rather unexpected, mechanism can be observed in the system depending on the choice of the material contrast ratio ρ. In particular, when ρ = 4, a dislocation pile up followed by dislocation reflection into Phase A is observed, as shown in Figs. 8b and 8d. Since both the QC fix as well as the FE-PN model were not set up for such a situation, i.e., no inclined glide planes or fully-resolved regions are present (recall Figs. 5 and 7b), they cannot account for dislocation reflection-unlike the fully flexible adaptive QC model.
It is worth noting that the presented speed-ups are obtained for two-dimensional atomic simulations only. In three-dimensional examples-which are more general and physically relevant-, a significantly larger number of atoms is present. The neighbor search algorithm within the considered cut-off radius will be significantly slower and the number of interatomic bonds is substantially higher. As a consequence, the QC simulations are expected to be much slower compared to the simulations of homogenized FE-PN model. Moreover, the presented example considers a simple Lennard-Jones pair potential. Use of a more advanced multibody potential, such as embedded atom model, will induce even further computing costs required by the QC simulations.
On the contrary, 3D extensions of the homogenized FE-PN model are expected to be relatively efficient provided that dislocation trajectories remain a priory known, and fixed glide planes can be considered. In a fully general case with arbitrary dislocation trajectories, however, sliding of glide planes should be considered in all possible lattice directions, which will result in small elastic regions representing individual atoms with multiple glide planes among them, each requiring its associated GSFE. Such a model may become computationally very expensive compared even to a full atomic simulation. In cases without any prior information on possible trajectories of dislocations, it is thus preferable to use adaptive QC approach instead.

Full lattice model
The second example considers a quasi-brittle (concrete) specimen subjected to a threepoint bending test, modeled at the mesoscale as a discrete lattice system shown in Fig. 12a. The domain Ω is of size L×H, inside which a homogeneous X-braced lattice with spacing l 0 is considered in between n par particles α positioned at r α stored in an index set N . All lattice interactions are modeled as damageable with an exponential softening law. To prevent spurious damage localization in the vicinity of prescribed displacements or applied loads, the lattice is made 100 times stiffer in padding regions with size 8l 0 × 6l 0 under the loading force and 22l 0 × 4l 0 around both supports. Upon loading, a localized crack growing along the symmetry plane from bottom to top is expected.
In analogy to the atomic system considered in the previous section, an inter-particle pair potential π αβ is introduced. It consists of two contributions, the elastically stored energy reflected by φ αβ , and a dissipation distance D αβ , Because dissipative mechanisms are included, a variational formulation of rate-independent systems is considered [47,48], and the interaction potential depends also on an internal variable ω αβ , reflecting the level of damage in the interaction connecting particles α and β. In addition, the pair potential is considered in an incremental fashion, i.e., at a time instance t k , π αβ k = π αβ (t k ), and depends on the configuration of the system in the current as well as the previous time step. To allow damage processes to evolve only under tension and not under compression, the elastically stored energy is split into two parts, (1−ω αβ )φ αβ (r αβ + )+φ αβ (r αβ − ), in which only the first term is affected by the damage variable. The two distance quantities, r αβ + = max (r αβ , r αβ 0 ) and r αβ − = min (r αβ , r αβ 0 ) (assuming φ αβ (r αβ 0 ) = 0), ensure that the damage variable weakens the interaction only under tension, i.e., for r αβ > r αβ 0 . A quadratic elastic potential is considered, expressed as a function of strain where, in analogy to the atomic system, r αβ = r β − r α is a vector of relative positioning, r αβ = r αβ 2 the Euclidean distance between a pair of particles α and β, and EA is the normal cross-sectional stiffness of the bond.
The dissipation distance D αβ measures the energy dissipated by a single interaction during the evolution of the damage variable between two consecutive states, ω αβ 1 and ω αβ 2 , i.e., where D(ω) is the energy dissipated during a unidirectional damage process up to the damage level ω. It is defined implicitly such that the following exponential damage law results where ε 0 is the limit elastic strain for which damage starts to evolve, and ε f characterizes the slope of the softening branch in the associated stress-strain diagram, and where the upper index αβ has been dropped for brevity; for more details see [62,Section 4.1]. The adopted constitutive parameters are specified in Tab. 5.   In analogy to Eq. (9), with the introduction of a time discretization for a considered time horizon t ∈ [0, T ] from Eq. (8), the evolution of the fully-resolved damageable lattice system is governed by the minimization of the total incremental energy Π k , minimized over all admissible configurations v with an initial condition q T 0 = [r T 0 , 0 T ], where q T = [r T , z T ] is a state variable of the entire system, consisting of the kinematic variable r storing the positions of all particles, whereas z T = [ω 1 , . . . , ω n int ] collects the damage variables of all n int interactions, and q k = q(t k ). The incremental energy of the entire system is obtained by collecting contributions from all interactions, i.e., where is the internally stored energy 2 , and is the global dissipation distance (recall also Eqs. (25) and (28)), B α is the initial set of nearest-neighbors associated with a particle α (cf. Fig. 12b), whereas f k is a vector of external forces applied at a time instance t k . At each time step, the minimization problem (30) is solved to obtain a local minimum. Because multiple minima may exist, the energetic solution corresponds to the one that satisfies the following energy balance which equates a sum of the internally stored energy and the dissipated energy, defined as with the work performed by the external forces To ensure that obtained minimum corresponds to a physically relevant energetic solution, the energy balance constraint of Eq. (34) is closely monitored in each time increment. Upon its violation, the considered time step is restarted with a smaller loading increment. To avoid such situations from occurring, the evolution path of the system is controlled using a dissipation driven arc-length method or crack mouth opening displacement control, as elaborated in more detail below in Section 3.5. For further details on the variational formulation of rate independent systems see, e.g., [49] (Section 3.2), and for applications to lattice systems, e.g., [39,60,27,61,62].

Quasicontinuum model
In analogy to Section 2.2, the two QC steps are introduced as follows. First, the interpolation step is introduced for the kinematic variable, cf. Eq. (11), while the internal variable z is not reduced for simplicity and to avoid issues with nonuniqueness in the prolongation operation and mesh coarsening as discussed in [60, Section 3.1]. The reduced QC state variable then reads q T red = [g T , z T ]. To allow for crack paths to be independent of the mesh topology, an extended variant of a QC methodology introduced in [61] is adopted. This extension allows for an efficient mesh coarsening through Heaviside step function enrichments used in the interpolation matrix Φ. The generalized kinematic variable g then stores the positions of all representative particles (r rep , analogous to atomic QC) as well as coefficients of linear combinations of the special function enrichments.
The summation step again samples the total incremental energy based on only a small set of selected sampling atoms stored in an index set N S , i.e., where w α and N S correspond to the central summation rule of Beex et al. [9], generalized for Heaviside type of enrichments according to [61]. The energy associated with a lattice site α is defined as π α k (q k ; q k−1 ) = 1 2 β∈Bα π αβ k (r αβ k , ω αβ k ; ω αβ k−1 ), (39) where the incremental energy of a single interaction π αβ k has been defined in Eq. (25). A QC method therefore minimizes the approximate total incremental energy of Eq. (38) with respect to the reduced state variable q red , mapped on the state of the entire system through the interpolation relation of Eq. (37).
The area of high interest, i.e., the fully-resolved region, is allowed to evolve adaptively throughout the simulation. Available options for the refinement criteria have been recently discussed in [22]. From the presented options, only a few studies presented concepts suitable for applications to structural lattices [46,62,61,51,22]. In this work, the following mesh error indicator, employed also in [62], is adopted for mesh refinement as well as for mesh coarsening. Interactions that are likely to be subjected to damage are identified through the refinement criterion, which compares the energy stored in an interaction (recall that damage is only allowed in tension) against a certain threshold energy, i.e., where θ r ∈ (0, 1) is a safety parameter and φ αβ th = φ αβ (r 0 (1+ε 0 )) is the elastic threshold energy at which the damage starts to evolve, and S K int denotes the set of all sampling interactions αβ located inside an element K. If at least one bond in S K int satisfies the condition of Eq. (40), element K is fully refined (i.e., all lattice sites inside the triangle are added as repnodes, and the triangulation is updated). Alternatively, element K is coarsened if the elastic energy of all bonds is lower than a given threshold, i.e., where θ c ∈ (0, 1) is a coarsening parameter satisfying θ c < θ r . Note that in Eq. (41) all bonds, i.e., including those in compression, are considered. The fully-resolved region is adaptively updated following the same three steps (i)-(iii) presented at the end of Section 2.2. No initial mesh refinement is needed, unlike in the atomic QC system, since no initiation analogous to the Volterra perturbation is employed (which needs to be resolved by the adopted QC system).

Finite element cohesive zone model
The initial lattice problem of Fig. 12a is translated into a continuous Finite Element Cohesive Zone (FE-CZ) model shown in Fig. 13, in which it is assumed that only one potential crack path is present in an otherwise continuous specimen. In contrast to a discrete lattice, the continuous domain does not require local stiffening under the applied boundary conditions to prevent excessive deformation. A cohesive zone Γ cz , located along the vertical axis of symmetry, splits the homogeneous domain Ω into two parts Ω L and Ω R , i.e., where x cz = 0 denotes the horizontal position of the cohesive zone. The mechanical behavior of the FE-CZ system is governed by the minimization of the total potential energy defined as where the second term on the right-hand side represents a concentrated force F in Fig. 13, distributed through tractions t over a small portion of the specimen boundary Γ t . Both subdomains Ω L and Ω R are purely elastic isotropic and the deformation related to the crack propagation is localized within the cohesive zone Γ cz . The total internal energy Ψ (cf. also Eq. (16)) reads where ψ e is the elastic strain energy density defined in Eq. (17), and ψ cz is the cohesive zone potential (defined solely on Γ cz ). Both are obtained through the homogenization of the underlying lattice as follows.

Calibration
For the calibration of the FE-CZ model, a single lattice unit cell as shown in Fig. 14a, is considered. The normal cross-sectional stiffness of the horizontal and vertical bonds are reduced to half of its true value to respect the considered periodicity. The corresponding assembled effective stiffness tensor D lat , listed in Tab. 6, is compared to an isotropic plane strain stiffness tensor D(E iso , ν iso ). The resulting effective elastic parameters are identified as where l z is the unit out-of-plane thickness. The fixed Poisson's ratio is given by the fact that D lat 1122 = D lat 1212 and the effective Young's modulus is obtained by fitting D lat 1111 = D 1111 (E iso , ν iso ). Note that for the considered setting of lattice structure (identical normal stiffness of all bonds), the resulting stiffness tensor D lat is not fully isotropic, cf., e.g., [50,Appendix B] for more details. The cohesive zone potential ψ cz is a function of the displacement jump across the cohesive zone Γ cz , see also Eq. (22). Due to the symmetry of the considered problem, the tangential component ∆ t vanishes, i.e., u L · e y = u R · e y , on Γ cz , and hence the cohesive zone potential ψ cz becomes a function of the crack normal opening ∆ n only. The cohesive zone potential is obtained numerically from the response of the lattice unit cell subjected to a horizontal stretch of ∆ n /l 0 with the boundary conditions shown in Fig. 14a. The corresponding traction acting on the cohesive zone as a function of the normal crack opening is defined as and is plotted for the lattice constitutive parameters of Tab. 5 in Fig. 14b.

Results and comparison
Due to the present damage processes, arbitrarily small increments in externally applied boundary conditions might lead to an uncontrolled evolution of the crack, and hence also of the internal variable z in both the fully-resolved and QC system (i.e., a snap-back may occur), which in turn might result in the violation of the energy balance constraint of Eq. (34). To avoid such situations, the evolution path is closely controlled using a dissipation driven arclength method, see, e.g., [74,33,45] for more details. For the FE-CZ approach, the standard Crack Mouth Opening Displacement (CMOD) control method [36] is sufficient to reliably compute a mechanical response of the system.
The mechanical behavior of the fully-resolved mesoscopic system of Section 3.1 (referred to as Full) is adopted as the underlying reference solution. Two QC models, as described in Section 3.2, are then considered: (i) one with a coarse triangulation of maximum element size H/4 (referred to as QC coarse), and (ii) one with a fine triangulation of maximum element size H/16 (referred to as QC fine), see Figs. 15a and 15b. The FE-CZ model, detailed in Section 3.3, is considered for only one triangulation with a fixed mesh. A maximum element size H/16 adopted far from the cohesive zone, which is gradually refined down to H/256 = l 0 towards the assumed crack path Γ cz , cf. Fig. 15c. A mesh convergence study has been performed to verify that the adopted element size provides converged results. The initial triangulations of all reduced models are shown in Fig. 15, whereas the corresponding number of DOFs are listed in Tab. 7. All three effective models are compared against the fully-resolved simulation in terms of force-displacement curves, peak forces, crack lengths, and crack opening profiles.
With the increasing external load, individual bonds near the bottom center region start to damage until localization occurs, forming a full central crack propagating almost vertically across the specimen's height, cf. Fig. 16b. The corresponding evolutions of the normalized magnitude of the externally applied force F/(EA) are shown in Fig. 17 as a function of the normalized vertical displacement measured under the force, u/l 0 . The specimen exhibits, after the initial linear response, a significant softening followed by a severe snap-back once the central crack localizes. While both QC approaches manage to capture the initial ductile behavior (the QC fine being almost indistinguishable from the Full solution), the FE-CZ completely omits this mechanism. Upon crack localization, however, all effective methods provide satisfactory qualitative description of the force-displacement curve. The QC coarse method suffers from an initially too stiff response as a consequence of a rather coarse mesh, although with adaptive refinement this discrepancy drops rapidly. To compare all methods quantitatively, the relative errors in the initial elastic stiffnesses and peak forces are summarized in Tab. 7. Here we conclude that the FE-CZ model underestimates the peak force by more than 23%, while being very accurate in the initial elastic stiffness. Although the QC coarse method overestimates the initial elastic stiffness by approximately 21%, it delivers an accurate peak force. The QC fine method provides very accurate results in both considered quantities.
The rather unexpected pronounced softening before crack localization is caused by an extensive damage region in which individual horizontal bonds are damaging or unload elastically in an alternating way, forming a checker-board pattern, cf. Fig. 16a. This region is hereafter referred to as the Distributed Damage Zone (DDZ). This partially damaged region expands progressively along the bottom edge of the specimen and upwards, featuring multiple possible crack paths (emerging in Fig. 16a) until localization occurs and the central crack is formed. A similar behavior is observed for both Full and QC models, whereas the FE-CZ model is unable to capture the initial DDZ as a consequence of the linear-elastic constitutive law used in the bulk material, which results in significant inaccuracies before the crack localizes (recall Fig. 17). Additional simulations (not shown) indicate that the observed errors further increase with increasing ductility of the considered exponential softening law of the individual bonds.
To verify that the distributed damage zone is at the root of the inaccuracy observed in the FE-CZ model, an additional simulation has been performed, in which damage was confined to evolve only in bonds located within 2 l 0 distance from the vertical axis of the symmetry, whereas the remainder of the specimen was considered as purely elastic. This version of the full model provided results (referred to as Lattice without DDZ) closely resembling those of the FE-CZ model. In particular, the relative error in the peak force of the FE-CZ model drops to 4.41% compared to this reduced version of the full model, thus revealing that the error in the peak force indeed results from the neglected distributed damage zone. To increase the accuracy of the FE-CZ model, multiple cohesive cracks should be realized in the distributed damage zone or a more involved constitutive model with damage should be employed for the bulk material.
To compare the crack lengths (Fig. 18) and the crack openings (Fig. 19), the central crack needs to be identified and its opening measured in discrete lattice systems. To this end, a fixed crack width of l 0 is considered, consisting of the most damaged bonds in each horizontal layer of bonds; if no damage occurs in a horizontal layer, the central bond (located at the specimen's vertical axis of symmetry) is used. The normal crack opening ∆ n is then defined as the difference between the horizontal displacements of the two end nodes of all crack bonds. Due to elasticity, negative crack openings are observed in compressive regions. The crack length, expressed as a function of the CMOD, is shown in Fig. 18. All approximate methods show comparable results, although both QC approaches are slightly more accurate, capturing also some of the initial irregularities of the Full model response. The FE-CZ method in general shows a smooth crack evolution as compared to the discrete models, initially overestimating the crack length, but then converging rapidly to the reference.
The corresponding crack opening diagrams, represented as the normalized normal opening ∆ n /l 0 and plotted against the normalized vertical coordinate y/l 0 , are shown in Fig. 19. The presented results correspond to a pre-peak configuration of F/EA = 0.24 (Fig. 19a), and to a post-peak configuration of F/EA = 0.035 (Fig. 19b). Here we clearly see that the FE-CZ and Full model without distributed damage zone (Lattice without distributed damage zone) tend to show significantly larger crack mouth openings before the crack localizes, whereas both QC models provide accurate predictions. Upon crack localization, in accordance with previously discussed results, all effective methods achieve a good accuracy.   The numerical performances of all methods are shown in Fig. 20 in terms of the number of DOFs relative to the number of DOFs of the fully-resolved system. Here we see that although both QC models start with noticeably different meshes (Figs. 15a and 15b), once the fullyresolved DDZ develops, the relative number of DOFs increases from approximately 0.5% to roughly 15% for both QC approaches. Upon crack localization, the mesh coarsens and the relative number of DOFs drops to 0.5% again. The FE-CZ model has no adaptivity, having a constant 1% relative number of DOFs. The achieved performance in terms of computing times is summarized in Tab 7. Here we see that in spite of relatively large savings in the number of DOFs, both QC methods attain a speed-up of only a factor of 2.5. This is caused mainly by the initial full refinement of the entire distributed damage zone (covering almost 1/7 of the domain in both QC models) and its subsequent coarsening, which require a substantial number of mesh iterations. This behavior can be avoided by using more aggressive mesh refinement/coarsening strategies, in which the mesh is updated only in selected steps or only a limited number of mesh iterations is allowed in each time step. Such modifications can significantly reduce the number of mesh iterations and thus speed up the simulation. However, if too aggressive, they may compromise the accuracy of the QC simulations. In cases with more localized phenomena, a higher computational gain by the QC methods is expected. The FE-CZ model, on the contrary, provides a substantial speed-up of the order of 50, but at the cost of a lower accuracy as compared to the Full solution.
Unlike the atomistic problem of Section 2, extensions to more relevant 3D structural lattices can be realized without significant changes to both the general QC as well as the homogenized model with fixed cohesive zones. In situations in which a crack trajectory is not known a priori and fixed cohesive zones cannot be used, alternative options for the homog-  enized model may be used. For instance, either multiple cohesive zones can be considered in the region in which crack is expected, or cohesive zones can be built in between all elements [78]. However, the latter option significantly increases the induced computing efforts and the resulting crack trajectory may moreover be mesh dependent. The extended finite element method, with cohesive formulation, is another interesting alternative [76,14,78], suffering, nevertheless, from technical difficulties for situations involving crack curving and crack branching.

Summary and conclusion
This paper has provided a detailed description and thorough comparison of two classes of homogenization techniques towards effective representation of two-dimensional discrete systems. The first class, referred to as QuasiContinuum (QC) based methods, considers the fully-resolved underlying discrete system, which is subsequently reduced through suitable mathematical tools such as projection and reduced integration. The second class, referred to as homogenization based methods, consists in homogenizing the underlying discrete system first, into which localized discrete mechanisms are subsequently embedded. The necessary theoretical basis for both classes of methods has been reviewed, and their performance demonstrated on two representative examples considered at the nano-and meso-scale, revealing their strengths and weaknesses.
The first class, i.e., QC based models, is capable of capturing all important phenomena of the underlying full model. The obtained results are thus sufficiently accurate for the examples considered at both scales, reproducing also unexpected mechanical behavior such as dislocation reflection and distributed cracking. At the same time, QC based models provide a significant reduction of the original problem in terms of the number of Degrees Of Freedom (DOFs), although the achieved savings in computing times might be less significant due to the computational cost involved with mesh adaptivity. Typical speed-ups achieved in this work are of the order of 2.5-15.
The second class, i.e., homogenization based approaches, represented by the Peierls-Nabarro model at the nanoscale and by the cohesive zone model at the mesoscale, is often more efficient as compared to the QC based methods in terms of computational speedup, although it is usually less accurate because the effects that have their origin in the discreteness of a system (e.g., Peierls stress) cannot be captured by the continuum approach. Typical savings in computing times corresponded to 5-50 in comparison with the full model. Although primal kinematic quantities such as overall displacement, dislocation positions, or crack opening and crack length are usually captured with adequate accuracy, the associated conjugate quantities such as transmission stress or maximum peak force suffer from more significant errors. The most significant weakness of the homogenization-based models is that they are constructed based on prior assumptions on the response of the system, thus being unable to capture general behavior and unexpected phenomena such as dislocation reflection, distributed cracking, or crack branching. This may be acceptable if one is certain that no such phenomena will occur, or even desired if one particular phenomena is to be studied that would only be clouded by the unexpected phenomena. If this is not the case, model limitations can be lifted by adopting, e.g., multiple glide planes of different orientations, assuming all inter-element interfaces as cohesive zones, or using an extended formulation combined with cohesive zone models. Such extensions might, nevertheless, significantly complicate the entire procedure and increase the associated computing time.
The resulting computational savings achieved by all reduced models are substantially affected by the level of scale separation, which is not very large in the presented representative examples. With increasing scale separation, e.g., increasing the size of specimen domain while fixing the lattice spacing, a significant improvement in observed speed-ups may be reached.
Both presented classes are sufficiently general tools for predicting localized phenomena in discrete systems. The homogenization based models are usually more efficient, although they might lack a sufficient level of detail and accuracy. The QC based methods, on the other hand, exhibit greater flexibility compensated by higher computing costs. Which method to adopt thus strongly depends on the desired accuracy and the acceptable computational costs. However, one should keep in mind the potential occurrence of unexpected phenomena. The QC based methods can thus be considered a safe option that may bring, nevertheless, only limited gain. Larger gain may be achieved by dedicated homogenized models, but safe option would be to check their results by a flexible QC method to verify occurrence of any unexpected phenomena.