A comparative accuracy and convergence study of eigenerosion and phase-field models of fracture

We compare the accuracy, convergence rate and computational cost of eigenerosion (EE) and phase-field (PF) methods. For purposes of comparison, we specifically consider the standard test case of a center-crack panel loaded in biaxial tension and assess the convergence of the energy error as the length scale parameter and mesh size tend to zero simultaneously. The panel is discretized by means of a regular mesh consisting of standard bilinear or Q1 elements. The exact stresses from the known analytical linear elastic solution are applied to the boundary. All element integrals over the interior and the boundary of the domain are evaluated exactly using the symbolic computation program Mathematica. When the EE inelastic energy is enhanced by means of Richardson extrapolation, EE is found to converge at twice the rate of PF and to exhibit much better accuracy. In addition, EE affords a one-order-of-magnitude computational speed-up over PF.


Introduction
The tracking of crack growth in solids is a free-discontinuity problem involving the formation of new internal surfaces. Modelling crack propagation remains a challenging problem in computational mechanics. In recent years, the phase-field (PF) method has gained in popularity, especially for crack propagation and tracking problems, cf., e. g., [1,2,3,4,5,6,7,8]. Proposed originally by Ambrosio and Tortorelli [9,10] as a regularization of the Mumford-Shah image segmentation functional, the PF method introduces an auxiliary continuous field, the phase field, as a means of representing the state of the material in vicinity of the crack. In effect, the phase field smooths the sharp surfaces of the cracks over a neighboring volume of finite thickness. In this way, the problem is reformulated in terms of displacement and phase fields defined over the entire volume of the domain. The governing equations define a system of second-order partial differential equations, thus in principle eschewing the difficulties inherent to evolving boundaries and discontinuities. However, the convenience of the initial implementation comes at the price of exceedingly fine discretization requirements in the vicinity of the cracks, a doubling of degrees of freedom, a sensitive dependence on the choice of intrinsic length scale, strongly non-linear and possibly unstable dynamics, difficulties enforcing strict irreversibility, no-healing and positive dissipation, difficulties enforcing crack closure, difficulties enforcing modemixity dependent fracture criteria, onerous computing time requirements and other difficulties, which need to be carefully addressed and assessed.
Element erosion (ER) methods [11,12,13,14], consisting of approximating cracks as notches of small but finite width, supply another wellestablished class of computational methods for simulating crack growth which has been extensively used to simulate fracture in a number of areas of application, including fragmentation and terminal ballistics. In seminal work, Negri [15] noted that some of the early versions of element erosion fail to converge, or converge to the wrong limit, due to mesh-dependency of the crack path, and provided a remedy based on the use of local averages over intermediate scales [16,17]. Subsequent enhancements of element erosion incorporating such local averaging [18,19,20,21,22,23,24,25,26] are provably convergent to Griffith fracture in the limit of vanishingly small mesh sizes [18].
Phase-field and element erosion methods have a common variational structure: i) an elastic energy-release mechanism, namely, progressive damage in the case of PF and abrupt damage in the case of ER; and ii) an energy cost of damage, derived from the phase-field and its gradients in the case of PF and from an estimate of the fracture area in the case of ER. In both cases, the static equilibrium configurations of the solid follow from global energy minimization. In addition, crack propagation is modeled in both cases by means of a rate-independent gradient flow that balances elastic energy-release rate and dissipation.
We begin by formalizing the commonalities between PF and ER methods and showing how they are special cases of a common variational structure based on the general notion of eigendeformation. Eigendeformations are widely used in mechanics to describe deformation modes that cost no local energy, cf., e. g., [27]. In some sense, eigendeformations provide the most general representation of elastic energy-release mechanisms. Not surprisingly, therefore, the method of eigendeformations provides a common framework for both PF and ER methods. Within this unified view, PF and ER methods simply correspond to particular choices of restricted classes of allowable eigendeformations and cost functions thereof.
Whereas the convergence properties of finite-element approximations of EE and PF is well established mathematically [28,18], a direct quantitative comparison and assessment of both methods appears to have been missing. In this work, we endeavor to fill that gap by means of selected numerical tests. By convergence we specifically understand convergence of the EE and PF solutions to the Griffith solution as the length parameter ǫ and the mesh size h both tend to zero. Thus, we regard the Griffith solution as exact and the EE and PF solutions as approximations thereof. Given the variational and energy minimization principles at work for both EE and PF, it is natural to measure errors and convergence rates in terms of energy and endeavor to ascertain the rate at which the EE and PF energies approach the limiting Griffith energy. We also recall that, for propagating cracks, the energy release rate supplies the requisite driving force for crack advance. Therefore, accuracy and convergence of the energy is a sine qua non prerequisite for the accuracy and convergence of the crack tracking problem.
It bears emphasis that we seek to characterize the convergence of solutions with respect to two parameters simultaneously, namely, ǫ and h. This double limit raises the fundamental question of the relative rates at which ǫ and h should be reduced to zero. Mathematical analysis [28,18] shows that convergence requires ǫ to decrease to zero more slowly than h, i. e., ǫ must be chosen on an scale intermediate between h and the size of the domain. Remarkably, this requirement is contrary to rules of thumb often used in practice that recommend setting ǫ to a fixed multiple of h. Here, we instead seek to optimize ǫ for given h as part of the approximation scheme. Specifically, for a given mesh size h we determine the optimal value ǫ h of ǫ by recourse to energy minimization, in the spirit of variational adaption. We show that ǫ h indeed yields the energy closest to the Griffith limit and thus minimizes the energy error for given mesh size h. The resulting meshsize convergence plots for EE and PF may therefore be viewed as the best possible for each method, which makes the method comparison fair.
We specifically consider the standard test case of a center-crack panel loaded in biaxial tension as a means of assessing the accuracy, convergence rate and computational cost performance of EE and PF. The panel is discretized by means of a regular square mesh consisting of standard bilinear or Q1 elements. In order to render the method comparison fair, all fields, including displacements and phase fields, are interpolated using the same shape functions. The exact stresses from the known analytical linear elastic solution are applied to the boundary. In order to deconvolve the discretization and quadrature errors, all element integrals over the interior and the boundary of the domain are computed exactly using the symbolic computation program Mathematica [29].
The results of the numerical tests reveal a superior accuracy and computational efficiency of EE over PF. In particular, when the inelastic EE energy is enhanced by means of Richardson extrapolation, EE converges at twice the rate of PF and exhibits better accuracy. In addition, EE is found to afford a one-order-of-magnitude computational speed-up over PF.

Multi-field models of brittle fracture
According to Griffith's criterion for fracture, in a brittle material crack growth is results from the competition between elastic energy minimization and the fracture energy cost of creating new surface. Assuming rate independence, crack growth in a solid occupying a domain Ω ⊂ R 3 is governed by the potential energy where is the total energy, including the elastic energy of the solid and the energy cost of fracture, W (ε(u)) denotes the strain energy density, ε(u) = sym ∇u the linearized strain tensor, u(x) the displacement field, dx the element of volume and the forcing terms (not spelled out for brevity) include body forces, boundary tractions and prescribed displacements. The jump set J u collects the cracks across which the displacement u may jump discontinuously and |J u | denotes the crack surface area. The material-specific parameter G c is the specific fracture energy density per unit area and measures the fracture strength of the solid.
The central and all-encompassing governing principle of energy minimization posits that the displacement field u at any given time is expected to minimize the potential energy Π(u) subject to monotonicity of the jump set J u , i. e., to the constraint that J u must contain all prior jump sets, and to crack closure constraints. In this manner, the problem of crack tracking is reduced to a pseudo-elastic problem, with monotonicity and closure constraints, for every state of loading. Such pseudo-elastic problems arise generally for rate-independent inelastic solids under monotonicity constraints (cf., e. g., [30] for a rigorous derivation) and were initially formulated in connection with deformation theory of plasticity [31].
The problem thus defined is a free-discontinuity problem in the sense that the displacement field u is allowed to be discontinuous and the discontinuity or jump set J u itself, i. e., the crack surface in the present application, is an unknown of the problem. The existence and approximation properties of such problems have been extensively investigated in the mathematical literature (cf., e. g., [32] for a review). Free-discontinuity problems are notoriously difficult to solve computationally, which has spurred the search for sundry regularizations of the problem that relax, to good computational advantage, the sharpness of the discontinuities. In the present work, we specifically focus on two such regularizations, eigenfracture and phase-field models, which we briefly summarize next.

2.1.
Eigenfracture. The method of eigenfracture (EF) is an approximation scheme for generalized Griffith models based on the notion of eigendeformation [18]. The approximating energy functional is assumed to be of the form where ε * is the eigendeformation field that accounts for fracture, E e (u, ε * ) is the elastic energy, E i ǫ (ε * ) is the energy cost of the eigendeformation, or inelastic energy, and ǫ is a small length parameter. The elastic energy E e (u, ε * ) follows as the integral over the entire domain of the strain energy density W as a function of the total strain ε(u) reduced by the eigenstrain ε * . In this manner, eigendeformations allow the displacement field to develop jumps at no cost in elastic energy.
This local relaxation comes at the expense of a certain amount of fracture energy. The challenge in regularized models of fracture is to estimate the inelastic fracture energy E i ǫ (ε * ) in a manner that converges properly as ǫ → 0. In the method of eigenfracture, the crack area is estimated as the volume of the ǫ-neighborhood {ε * = 0} ǫ of the support {ε * = 0} of the eigendeformations scaled by 1/ǫ, cf. Fig. 3b. Specifically, in this construction {ε * = 0} is the set of points where the eigendeformations differ from zero, {ε * = 0} ǫ is the ǫ-neighborhood of {ε * = 0}, i. e., the set of points at a distance to {ε * = 0} less or equal to ǫ, and |{ε * = 0} ǫ | is the volume of {ε * = 0} ǫ .
Remarkably, the method of eigenfracture is provably convergent [18], in the sense that the total energy (3) Γ-converges to the Griffith energy (2) in the limit of ǫ → 0. This convergence property shows that the eigenfracture method is indeed physically and mathematically sound. We recall that Γconvergence of the energy functionals in turn implies convergence of the solutions as ǫ → 0, i. e., the eigenfracture solutions converge to the solutions of Griffith fracture in the limit of vanishingly small length parameter ǫ.

2.2.
Phase-field models of fracture. In the PF approximation of Griffith fracture, the state of the material is characterized by an additional continuous field v(x) taking values in the interval [0, 1] and v = 0 at the crack. The crack set J u is then approximated as a diffuse interface where v = 1. The corresponding fracture model traces back to the pioneering work of Ambrosio and Tortorelli [10], who showed that a two-field functional Γ-converges to the Mumford-Shah functional of image segmentation. Generalized to threedimensional elasticity, the two-field functional of Ambrosio and Tortorelli assumes the form where ǫ is a small length parameter and o(ǫ) stands in for a positive function that decreases to zero faster than the small parameter ǫ. The work of Ambrosio and Tortorelli, and other similar works [33,32], subsequently spawned numerous variants, extensions and implementations (e. g., [34,35,1,15,16,36], but the differential structure of the fracture energy E i ǫ (v) in (4) has remained essentially unchanged in the later works.

2.3.
Eigenerosion. Eigenerosion (EE) supplies an efficient implementation of the eigenfracture model [19]. To establish the connection between eigenfracture, eq. (3), and eigenerosion, assume that W (ε) is quadratic and restrict eigendeformations to the particular form with w taking the values 0 or 1, i. e., w(x) ∈ {0, 1}. Inserted into eq. (3) this gives the EE functional By Jensen's inequality and properties of extreme points [37], it follows that the range of w can be extended to the entire interval [0, 1], i. e., 0 ≤ w(x) ≤ 1, without changing the solutions. It thus follows that EE is a restricted form of eigenfracture and, therefore, it supplies an upper bound of the eigenfracture energy in general. Evidently, the EE energy (6) may be regarded as a PF model with phase field and a fracture energy computed by the ǫ-neighborhood construction. Conversely, PF models of fracture may be viewed as special cases of EE, and hence eigenfracture, where the fracture energy is of the Ambrosio-Tortorelli type.
The great advantage of the EE model (6) vs. the conventional Ambrosio-Tortorelli-type phase-field model (4) is that in the former, eq. (6), the phasefield is undifferentiated and evaluates the fracture energy through an integral expression, whereas the latter, eq. (4), requires the phase-field to be differentiated. Differentiation in turn requires regularity and conforming interpolation, e. g., by the finite-element method. By contrast, the integral form of the fracture energy in (6) allows the phase-field to be approximated, e. g., as piecewise constant 0 or 1, which leads to a considerable increase in implementational simplicity and robustness [19,20].
2.4. Non-local fracture as an Artificial Neural Network. Artificial neural networks provide a compelling interpretation of the ǫ-neighborhood construction that suggests an entire class of extensions thereof. To make this connection, we introduce the mollifier and the activation function Next we note that the function obtained by taking the convolution of ϕ ǫ and w, is positive in the ǫ-neighborhood {w = 0} ǫ of the crack set and vanishes elsewhere. Therefore, the filtered function f (w ǫ (x)) is 1 in the ǫ-neighborhood {w = 0} ǫ and vanishes elsewhere, i. e. it is the characteristic function of the ǫ-neighborhood. Finally, we have which supplies an integral representation of the ǫ-neighborhood construction.
In (11), we immediately recognize the structure characteristic of an artificial neural network (cf., e. g., [38]), Fig. 1. Thus, the fracture energy E i ǫ (w), which is the output of the network, follows from the input w through the composition of three operations. The first operation is a convolution ϕ ǫ * w, which defines a linear neural network. The outcome w ǫ of this operation is filtered locally by means of the activation function f in the form of a binary switch, with the result that f (w ǫ ) is the characteristic function of the ǫ-neighborhood of the crack set. Finally, the fracture energy follows as an integral of f (w ǫ ), suitably scaled by G c /2ǫ.
The artificial neural network interpretation suggests an extension of eigenfracture to a more general class of fracture models where the fracture energy is of the form (11) but ϕ ǫ and f and a general mollifier and activation function. This generalization immediately raises the question of how the accuracy and convergence properties of the model depend on the choice of ϕ ǫ and f .

Test case: Slit crack under all-around tension
We proceed to assess the accuracy and convergence of PF and EE approaches by recourse to the standard test case of a slit crack in an infinite solid subjected to all around tension, Fig. 2. We regard the Griffith solution as exact and the EE and PF solutions as approximations thereof. Since both the PF and EE problems are energy minimization problems, we specifically monitor energy errors and seek to characterize the convergence with respect to the length parameter ǫ and mesh size h simultaneously. 3.1. Exact reference results. We specifically consider an infinite solid containing a straight crack of length 2a deforming in plane strain under the action of equibiaxial stress σ 0 at infinity, see Fig. 2(a). The crack-tip stress field is mode I since the loads are symmetric with respect to the crack line.
Conveniently, the problem has an exact solution [39], namely, where the coordinates θ, θ 1 , θ 2 , r, r 1 , and r 2 are defined in Fig. 2(b) and the lower indices correspond to cartesian coordinates (x 1 , x 2 ) aligned and centered at the crack. Assuming plane strain conditions, the strains follow from Hooke's law as where E denotes the Young modulus and ν the Poisson's ratio. The displacements u can then be computed by integrating the relations (14) ε 11 = ∂u 1 ∂x 1 , using Cesaro's method. Finally, we recall that the strain-energy density of the solid is (15) .
In the absence of the crack, i. e., for crack length a = 0, an application of Clapeyron's theorem gives (20) Π where |Ω| is the area of Ω and W 0 is given by (18). For a finite crack, the minimum value of the potential energy follows directly from an application of Rice's J-integral [40]. Specifically, the energy-release rate is given by where Γ denotes a counter-clockwise closed contour surrounding the crack and contained in Ω, n is its outward unit normal, and ds is the element of arc-length on Γ. Choosing Γ to coincide with the flanks of the crack, together with small loops at the tips, and using the asymptotic K-field gives where K I is the mode I stress-intensity factor, which can be computed directly from the stress field (12), with the result Inserting (23) into (22), integrating with respect to a, and using (20) with (18) as initial condition gives In addition, according to the Griffith model (2), the inelastic or fracture energy expended in the extension of the crack is The exact results (24) and (25) are subsequently taken as a convenient basis for the analysis of the accuracy and convergence of EE and PF approximation schemes.

Discretization and implementation
Next we describe the discretization of the EE and PF models used in calculations. All calculations are performed on a square domain Ω of size D ≫ 2a.
with the energy E ǫ (u, w) as in (6), is discretized by means of a regular mesh consisting of standard bilinear or Q1 elements of size h. The exact stresses σ ij (x) in (26) are taken directly from the analytical solution (12). In order to disentangle the convergence with respect to the mesh size h from quadrature error, in all the calculations we evaluate all element integrals over the interior and the boundary of the domain exactly using the symbolic computation program Mathematica [29]. The implementation of the EE fracture energy is illustrated in Fig. 3. The crack is represented as a string of missing, or 'eroded', elements where w = 0, dark gray elements in Fig. 3(a), with w = 1 elsewhere. The eroded elements approximate the crack geometry and the corresponding inelastic or fracture energy is approximated as where ǫ is a length parameter and A ǫ,h is the area of the ǫ-neighborhood of the eroded elements, light gray area in Fig. 3(b), i. e., the set of points at a distance smaller or equal to ǫ from the eroded elements. Examples concerned with crack growth through slanted meshes have been presented in [19]. For a slit crack centered and aligned with the mesh, the eroded elements are precisely those which contain the crack, and A ǫ,h follows exactly as (28) A ǫ,h = ⌈2a/h⌉h 2 + 2(⌈2a/h⌉ + 1)hǫ + πǫ 2 .
As expected, the EE approximation of the inelastic energy E i EE,ǫ,h depends on the choice of length parameter ǫ. From a variational viewpoint, the optimal value ǫ h of ǫ is that which minimizes E i EE,ǫ,h , namely, The corresponding optimal inelastic energy is an asymptotic expansion of (30) gives in the limit of h/a ≪ 1 A similar asymptotic expansion of the optimal inelastic energy (31) likewise gives We observe from (33) that, to leading order, the optimal length parameter ǫ h scales as √ 2ah, i. e., as the geometrical mean of the crack length and the mesh size. Thus, the optimal length parameter depends not only on the mesh size but also on the geometry of the crack and, in general, of the domain. We note that ǫ h → 0 as h → 0, albeit at a slower rate, as required by convergence [18].
We also observe from (34) that, to leading order, the fracture energy error is of order O(h 1/2 ). This rate of convergence is slower than the O(h) rate of convergence of the elastic energy and, therefore, dominates the overall energy error. However, this loss of convergence can be remedied simply by recourse to a standard Richardson extrapolation technique (cf., e. g., [41,42]). To this end, we note from (34) that the fracture energy attendant to a mesh of size 2h is We may replace E i EE,h by the weighted sum without disturbing the limit. We then choose the weight λ so as to cancel the O(h 1/2 ) term, with the result From (34) and (36) it then follows that , as desired. Inserting (31) and (37) into (36), we find explicitly. This simple Richardson extrapolation effectively eliminates the low-order accuracy of the original ǫ-neighborhood construction and restores the full order of convergence expected of the finite element method.

4.2.
Phase-field model. In order to have a fair comparison with EE, we discretize the potential energy (4), by means of a regular mesh likewise consisting of standard bilinear or Q1 elements of size h. Conforming interpolation is used for both the displacement and the phase fields. In order to separate interpolation errors from numerical quadrature errors, we again evaluate all element integrals over the interior and boundary of the domain exactly using the symbolic computation program Mathematica [29]. The phase field is unconstrained and, therefore, satisfies free Neumann boundary conditions on the outer boundary of the domain. The unknown fields (u h , v h ) are solved iteratively by the method of alternating directions [43], i. e., by successively fixing one of the fields and solving for the other. Conveniently, the scheme reduces the solution to a sequence of linear problems (cf., e. g., [44,45]). The iteration is primed by setting, as initial condition for the iteration, v h = 0 on the nodes lying on the crack and v h = 1 elsewhere. The iteration may then be expected to approximate the solution for a crack of length 2a provided that the applied stress σ 0 equals the critical stress for crack extension, i. e., if as in this case the exact elasticity solution (12) and the crack length 2a jointly minimize the Griffith potential energy (1).
The choice of the length parameter ǫ, and its relation to the mesh size h and geometrical features of the domain, is known to have a strong effect on the accuracy and convergence properties of PF approximations [28]. In particular, convergence requires h to decrease to zero faster than ǫ ( [28], Theorems 4.1 and 5.1). We expect the exact potential energy (24) to be approached by the PF potential energy from above (cf., e. g., Fig. 4b). For fixed h, we also expect the PF potential energy to exhibit a minimum at a certain value ǫ h of ǫ. Thus, for small ǫ the mesh size h is unable to resolve the width of the crack, resulting in an overly stiff response and high PF potential energy. Contrariwise, since the exact potential energy (24) is attained from above for ǫ → 0, the PF potential energy diverges from the exact value for large ǫ. As a consequence of these opposing trends, for fixed h the PF potential energy attains a minimum at a certain ǫ h , as surmised (cf., e. g., Fig. 4b). Since, as already noted, the exact potential energy (24) is approached by the PF potential energy from above, the energy minimizing ǫ h results in the least energy error for given h and is therefore optimal from the standpoint of convergence.
Evidently, ǫ h depends on the geometry of the crack and of the body, the state of damage and the discretization. In calculations, we determine ǫ h numerically by computing the PF minimum potential energy for given h over a range of equally-spaced values of ǫ, interpolating the computed energies as a function of ǫ and computing the minimum of the interpolated function (cf., e. g., Fig. 4b).
We proceed to investigate the relative accuracy and convergence rates of the EE and PF methods by way of numerical testing. To this end, we fix the domain size at D = 5, the crack length at 2a = 0.403215 and consider a sequence of five uniform meshes of sizes h/D = 0.02, 0.01, 0.005, 0.0025, and 0.00125. The numerical parameters used in calculations are listed in Table 5. We note that the parameters are chosen so as to satisfy the criticality condition (41).

5.1.
Optimal choice of length parameter. Fig. 4 shows the computed dependence of the minimum potential energy Π ǫ on ǫ at fixed h for both EE, eq. (26), and PF, eq. (40). As surmised, at fixed h the minimum potential energy Π ǫ attains a minimum at a well-defined optimal value ǫ h for both EE and PF. For EE, ǫ h is given analytically and in close form by (30). For PF, ǫ h is determined numerically by interpolating the results shown in Fig. 4b.  As expected [28,18], the optimal length parameter ǫ h converges to zero more slowly than the mesh size h. The scaling law (42) follows heuristically if we assume that near the crack the phase field varies on the scale of ǫ 2 /L, where L is a characteristic size (intrinsic geometric feature, e. g., domain size, crack size, ligament size). In order to resolve this variation, we must choose h ∼ ǫ 2 /L, whence (42) follows. The square-root nonlinearity of the scaling law (42), or equivalently ǫ h ∼ √ Lh, results in optimal values of the length parameter that may be strongly incommensurate with h, contrary to standard computational practice. Thus, for fine meshes, h ≪ L it follows that ǫ h ≫ h, i. e., ǫ h lags behind h and h resolves ǫ h finely. Contrariwise, for coarse meshes, h ≫ L, we have ǫ h ≪ h, i. e., ǫ h runs ahead of h and is unresolved by h. This latter regime is clearly visible in Fig. 4b, e. g., at h = 0.1, for which the corresponding ǫ h is over one order of magnitude smaller. 5.2. Energy convergence. Potential energy convergence plots for eigenerosion (EE), eigenerosion with Richardson extrapolation (EE+RE) and phase-field (PF) are shown in Fig. 6a. The plot displays least-square fits of the data by functions of the form (43) inf Π ǫ h = Π 0 + Ch α , with respect to the parameters C and α. The limiting potential energy Π 0 in (43) is given by (24). The exponent α is the rate of convergence and the constant C shifts the error line vertically in log-log convergence plots. All energy errors are computed using the optimal ǫ h corresponding to each mesh size, computed as described earlier. The figures visualize the least error in minimum potential energy as a function h for each method. All terms are normalized: the mesh size h is normalized by the crack length, 2a, and the potential energies are normalized by Π 0 . As may be seen from Fig. 6a, both PF and EE exhibit sublinear (squareroot) energy convergence, though the constant for EE is nearly one order of magnitude smaller (better) than that for PF, indicating superior accuracy of EE over PF under the conditions of the test. By far the best accuracy and convergence rate is obtained for EE+RE, which exhibits linear energy convergence, or twice the rate of convergence of EE and EF, and a better constant than that of both EE and PF. Fig. 6b shows a comparison between execution times for EE and PF. We recall that the preceding convergence calculations are carried out using exact integration of the element energy and nodal forces. This procedure has the advantage of singling out interpolation errors and deconvolving them from other sources of error such as numerical quadrature. However, the valuation of the exact integrals is costly and deviates from common practice, which invariable relies on numerical quadrature as a further approximation. Therefore, in order to obtain practical estimates of execution times, we repeat all calculations using standard finite element numerical quadrature. Specifically, we use four-point Gaussian quadrature for the element integrals and two-point Gaussian quadrature for the boundaryedge integrals. All calculations are performed using a sparse linear solver [46] in memory-shared configuration, using a single node, 16-core Intel Skylake (2.1 GHz), with 192 GB Memory and 2666 MT/s speed. Fig. 6b shows that, under the conditions of the test and for the particular computer architecture used in the calculations, EE is about one order of magnitude faster than PF. The higher efficiency of EE over PF is in fact expected, since EE entails displacement degrees of freedom and a no-iteration direct solution, whereas PF entails both displacement and phase-field degrees of freedom and an iterative solution.

Summary and concluding remarks
We have presented a comparison of the accuracy, convergence and computational cost performance of the eigenerosion (EE) and phase-field (PF) methods for brittle fracture. Both approaches operate on the principle of minimization of a potential energy functional that accounts for the elastic energy of the system, the inelastic energy attendant to crack growth and the work of the applied loads. Both approaches can be derived as special cases of the general method of eigendeformations by effecting particular choices of the eigendeformation field. The energy functionals incorporate an intrinsic length scale ǫ representing an effective crack width. The solution for Griffith fracture is obtained in the limit of ǫ → 0.
Whereas the convergence of finite-element approximations of EE and PF is well established mathematically [28,18], a head-on relative assessment of both methods appears to have been missing. For the standard test case of a center-crack panel loaded in biaxial tension, the results of the numerical tests reveal a superior accuracy and computational efficiency of EE over PF. In particular, when the accuracy of the EE inelastic energy is enhanced by means of Richardson extrapolation, EE+RE converges at twice the rate of both PF and the original low-order version of EE. In addition, EE affords a one-order-of-magnitude computational speed-up over PF.
There are other intangible benefits that confer EE an advantage over PF. Thus, element erosion is exceedingly easy to implement in accordance with a critical energy-release rate criterion and it guarantees automatically both irreversibility and positive dissipation. Another considerable disadvantage of PF relative to EE is that it requires iteration, a doubling of the degrees of freedom and defines an extremely nonlinear and non-convex problem with vastly many local minima, with the attendant difficulties in terms of stability and convergence.
In closing, we emphasize that the conclusions reached in this study are based on a limited and selected set of numerical tests. Additional tests would be greatly beneficial and are likely to shed further light on the relative merits of the methods.