Mesh-Based and Meshfree Reduced Order Phase-Field Models for Brittle Fracture: One Dimensional Problems

Modelling brittle fracture by a phase-field fracture formulation has now been widely accepted. However, the full-order phase-field fracture model implemented using finite elements results in a nonlinear coupled system for which simulations are very computationally demanding, particularly for parametrized problems when the randomness and uncertainty of material properties are considered. To tackle this issue, we present two reduced-order phase-field models for parametrized brittle fracture problems in this work. The first one is a mesh-based Proper Orthogonal Decomposition (POD) method. Both the Discrete Empirical Interpolation Method (DEIM) and the Matrix Discrete Empirical Interpolation Method ((M)DEIM) are adopted to approximate the nonlinear vectors and matrices. The second one is a meshfree Krigingmodel. For one-dimensional problems, served as proof-of-concept demonstrations, in which Young’s modulus and the fracture energy vary, the POD-based model can speed up the online computations eight-times, and for the Kriging model, the speed-up factor is 1100, albeit with a slightly lower accuracy. Another merit of the Kriging’s model is its non-intrusive nature, as one does not need to modify the full-order model code.


Introduction
Fracture is one of the most commonly-encountered failure modes of engineering materials and structures. Prevention of cracking-induced failure is, therefore, a major constraint in engineering designs. As with many other physical phenomena, computational modelling of fracture constitutes an indispensable tool not only to predict the failure of cracked structures, for which full-scale experiments are either too costly or even impracticable, but also to shed light onto understanding the fracture processes of many materials such as concrete, rock, ceramics, metals, biological soft tissues, etc. Within the context of continuum modelling of brittle fracture, this paper presents mesh-based and mesh-free reduced-order phase-field models. The proposed models can be used for parameter sensitivity analysis for brittle fracture problems, parameter estimation for phase-field brittle fracture vectors. The Kriging model, as a promising alternative approach, is a meshfree surrogate model that adopts statistical methods, which can provide deeper insight into the relationship between some outputs of interest and input design variables. Compared with POD-(M)DEIM, the advantage of the Kriging model is its fast online computations and lower computer storage.
In order to avoid all the complexities of PFMs so as to focus on the ROM itself, we have selected a one-dimensional (1D) PFM for quasi-static small strain brittle fracture. In this simple setting, one does not have to deal with strain energy splits, damage boundedness, irreversibility, and so on. We emphasize that the aims of this paper are not to solve the 1D parametrized brittle fracture problem by using a ROM, the solution of which can be found analytically [55], but to demonstrate how to build ROMs for PFMs and present preliminary 1D results. The proposed ROMs are inherently multi-dimensional. The fact that we have applied them to a very simple problem of a softening bar is due to the lack of computing resources to perform the intensive offline calculations (one might need to carry out 5000 2D/3D fracture simulations). However, we are not clear if ROMs can perform well for two-and three-dimensional PFMs.
Based on computations of a one-dimensional bar with varying the Young's modulus and fracture energy (thus, geometry, loadings, and boundary conditions are fixed even though they can also be parameters), the POD-(M)DEIM ROM can speed up the online computations eight-times, whereas for the Kriging model, the speed up factor is 1100, albeit with a slightly lower accuracy. Moreover, the Kriging model has the extra advantage of its non-intrusive nature in the sense that one does not need to modify the full-order model code. Needless to say, all these savings in the computational cost are achieved with extensive offline computations using the full model. These encouraging one-dimensional results are simply a proof-of-concept demonstration and serve as a platform to build ROMs for three-dimensional brittle fracture problems.
The remainder of this paper is structured as follows. Section 2 presents the formulation of the selected phase-field brittle fracture model, which includes the governing equations, the weak form, and the finite element solver. This is followed by Section 3, which is devoted to the presentation of the two reduced-order models: the POD-(M)DEIM ROM in Section 3.1 and the Kriging model in Section 3.2. Numerical examples are provided in Section 4 to assess the performance of these models. Conclusions and further works required to lift the limitations of the current work are given in Section 5. The POD algorithm is presented in Appendix A, and the analytical solution of the investigated model is given in Appendix B.

Phase-Field Model for Quasi-Brittle Fracture
This section briefly recalls the one-dimensional phase-field model for brittle fracture. Governing equations are given in Section 2.1, and the weak form and finite element discretization are presented in Section 2.2. We refer to [3,6,23] for details on various phase-field models for brittle and cohesive fracture.

Governing Equations
Let us consider a bar of length L, which is fixed at the left end (x = 0) and pulled at the right end (x = L). Without loss of generality, a unit cross-sectional area is assumed. The displacement and crack phase-field (or damage field) are represented by the functions u := u(x, d; E 0 , G c , b, u * ) and d := d(x, u; E 0 , G c , b), respectively. The following governing equations and boundary conditions hold [23]: u(x = 0) = 0, u(x = L) = u * , (essential boundary condition for u) d(x = 0, L) = 0, (essential boundary condition for d) where (x) and σ(x) are the strain and stress fields, respectively; E 0 and G c denote Young's modulus and the fracture energy of the material; b is the length scale introduced to approximate a sharp crack by a diffuse damage band; the loading is described via the imposed displacement u * . The essential boundary conditions for d are to prevent damage from initiating at either end of the bar. For this particular phase-field model, the damage boundedness condition 0 ≤ d ≤ 1 is automatically satisfied [6]. As only monolithic loadings are herein considered, damage irreversibilitẏ d ≥ 0 is also fulfilled without any special consideration. Note however that, if needed, it can be enforced quite straightforwardly using various techniques; see [23] for details. As can be seen, PFM involves the solution of two coupled partial differential equations: the equilibrium equation and the damage evolution equation. This usually leads to a misunderstanding that PFM is just another gradient-enhanced damage model developed by [56]. It has been computationally shown in [5,11] that when the length scale approaches zero, the PFM approaches a fracture mechanics model rather than a damage model. Note that body forces are omitted for simplicity. As can be seen from the non-homogeneous Dirichlet boundary condition u(x = L) = u * , we adopt a displacement control to trace the entire equilibrium path as snap-back does not occur for problems considered in this work. If needed, arc-length control can be used; see, e.g., [57].

Weak Forms and Finite Element Implementation
The weak form of the above governing equations is given by: for test functions δu and δd.
In order to simplify the notation, let us denote the vector of state variables by x = (u, d) and the vector of model parameters by µ = (L, E 0 , G c , b, u * ) ∈ Ω. This weak form is discretized by using standard finite elements, resulting in the following discrete equations (see Remark 1): where a(x; µ) ∈ R N andā(x; µ) ∈ R N are the nodal displacement vector and damage vector, respectively. K aa (x; µ) ∈ R N×N and K dd (x; µ) ∈ R N×N are the matrices, and f(x; µ) ∈ R N is the (external) force vector (we admit this terminology is not entirely correct, but as there is only one vector in the system, we hope that it does not cause any confuse). Here, N is the number of grid points. Note that one has to solve the displacement equation with the constraint a N = u * and the damage equation withā 1 =ā N = 0. The matrices and the force vector in the above are given by: where quantities with a bar indicate fixed values (as a staggered solver is used, which will be shortly discussed); N is the row vector of the shape functions; and ∇N denotes the row vector of the first derivatives of the shape functions. The symbol T denotes the transpose operator.

Remark 1.
Due to the non-convexity of the energy functional in terms of both displacements and damage, monolithic solvers are not very robust. That is why staggered solvers or alternating minimization solvers, where the off-diagonal coupling matrices are not needed, are popular in phase-field fracture [3,6].
The system of Equations (9)-(10) is solved using a staggered solver, also known as the Alternating Minimization (AM) solver. That is, one fixes the damage in the equilibrium equation and solves for the displacement. Next, the updated (latest) displacement field is used to calculate the driving force and substituted into the damage evolution equation to solve for the damage field. The process is repeated until convergence. These steps are summarized in Algorithm given in Box 1, where the notation a k n+1 signifies the nodal displacement vector at load increment n + 1 and AM iteration k; and a n+1 denotes the converged displacement vector at step n + 1. Basically, there are two computational bottlenecks: 1. the solution of two N × N systems for each AM iteration k and 2. the evaluation of the force vector f and matrices K aa and K dd .
which render PFM computations time-consuming and may not be feasible in situations where they have to be repeatedly executed a large number of times. Reduced-Order Models (ROM) can be applied for such problems to obtain a lower order efficient model. They are introduced subsequently in Section 3. Note that we do not focus on the cost of the robust-but-slow AM solver and refer to [23] for a discussion on efficient solvers for phase-field models. We refer to Remark 2 for extension to 2D/3D problems. Box 1. Quasi-static brittle PFM: AM solver for load step n + 1.

Reduced-Order Modelling
This section presents two reduced-order models, one based on the mesh-based POD method presented in Section 3.1 and the other based on a meshfree approach (cf. Section 3.2).

Mesh-Based Approach
Essentially, a ROM is carried out in two phases: a computationally-expensive offline phase and a computationally-efficient online phase. In the offline phase, a set of samples is collected from a standard analysis of the full-order model (in this context, the PFM simulation). This information is employed to construct a reduced-order model that is used in the online phase. In practice, the POD is applied to a set of samples collected from a full-order PFM analysis to compute a set of basis vectors (Section 3.1.1). These basis vectors are later used in the POD-DEIM to build the approximation for the force vector f (cf. Section 3.1.2) and in the POD-(M)DEIM for the matrices K aa and K dd (Section 3.1.3).

Parameterized and Nonlinear ROM Based on the Projection Framework
Consider the system of equations in Equations (9) and (10). An ROM of this system can be obtained by approximating the full-state vectors a andā as a linear combination of m andm basis vectors as follows, where a r ∈ R m andā r ∈ Rm are the reduced-order versions of the displacement and damage field, ∈ R N×m are the orthonormal bases. Now, at every iteration, step k in Box 1, projecting the system of Equations (9) and (10) onto the reduced spaces formed by these bases yields the lower order model as follows: where the reduced-order matrices and vector are given by: The ROM task is to find the bases V andV so that m N andm N, then to solve the system of Equations (15)-(19) using Box 1. How to determine these bases is subsequently discussed in Box 3. This task is simple and straightforward when the original system is an affine and linear system; and it would be complex for nonlinear and coupled systems. The construction and solution of the system (15)- (19) over previous solutions, i.e., (x k = (Va k r ,Vā k r )), and the parameter space are nontrivial. For instance, at each AM iteration k, the ROM requires firstly the re-construction of the full-order system matrices K aa , K dd and vector f corresponding to the parameters and previous solutions, which still depend on the dimension of the full model, and secondly, the projection of these matrices/vectors on the reduced spaces to obtain the corresponding reduced-order matrices and vector. This pure POD model may result in a longer and much more complicated computation than the original FOM. In this study, we implemented DEIM and (M)DEIM [31,35] in the offline stage to get rid of the full-dimension dependence of the ROM matrices and vector. Precisely, DEIM was used to approximate the nonlinear external force, that is the right-hand side of the damage sub-problem equation (Section 3.1.2), and (M)DEIM was adopted to approximate the nonlinear matrices (Section 3.1.3).

DEIM
The DEIM was applied to approximate the nonlinear function of the external force in Equations (19) and (13). The idea is to select sampling of the nonlinear terms combined with interpolation among these samples to recover an approximate nonlinear evaluation, as follows: where matrix Φ is the POD basis vectors of the nonlinear snapshots obtained from the FOM, . . , e n f ] ∈ R n×n f is the n f -indices matrix provided by DEIM, which we used here as the original proposed in [31], where e i = [0, . . . , 0, 1, 0, . . . , 0] T ∈ R n is the i th column of the identity matrix I n ∈ R n×n for i = 1, . . . , n f . Note that Φ is acting as a projector of the basis on vector f (Equation (21)), while P is acting as a filter matrix that returns the non-zero components of f (Equation (22)). In other words, it is used to define the reduced mesh (see Section 3.1.3). Let us define: The POD-DEIM reduced-order of the external force in (19) now has the form: The complexity of the evaluation of f r is now reduced to the evaluation of F r and a matrix multiplication (note that the computation of D was carried out in the offline phase). The POD-DEIM basis vectors were obtained using the algorithm in Box 2 where N s represents the number of sampling points (the number of points that discretize the parameter space), and N T , in this paper, denotes the number of load increments. We refer to Box A1 for the algorithm of POD(X f , f ). Note that the DEIM was used in the offline phase (see Box 3) to build the reduced matrices and vectors for the approximation of the FE external force and matrices.

Box 2. DEIM algorithm.
input : This max function returns the maximum value ρ at location Φ 1 . b e i is the i th column of an identity matrix.

(M)DEIM
The (M)DEIM was applied to approximate the nonlinear matrices in Equations (17) and (18). The idea is to express the nonlinear matrices (i.e., K dd (x k ; µ)) in vector format, then apply the DEIM to approximate it. Without loss of generality, let us consider K dd (x k ; µ) (the matrix of the damage sub-problem) and define the vector version of it: that is, k is obtained by stacking the columns of K dd . Then, this vector is approximated as: Here, Φm ∈ R N 2 ×mā is a nonlinear-parameter-independent basis and θm(x k ; µ) ∈ R mā is the coefficient vector.
Apply DEIM in Box 2 to a set of snapshots

obtain the basis Φm and the interpolation indices
The reduced vector is obtained by the following projection: Let us define Dm := (V T ⊗V T )Φm, which is a precomputed matrix; the online computation of the reduced vector becomes: where: Then, the (M)DEIM approximation matrix K dd (x k ; µ) is obtained by reversing the vec(·) operation. The crucial step in the online evaluation of k r is the computation of k I (x; µ). Within the framework of the finite element method, a reduced mesh concept, also called the reduced integration domain, can be used to evaluate k I (x; µ) efficiently. The basic idea is to loop over only elements belonging to I that contribute to the stiffness matrix. For more details of this technique, we refer to [35,59]. The exact same algorithm applies for K aa (x k ; µ) by defining k(x k ; µ) := vec(K aa (x k ; µ)) ∈ R N 2 and replacingm by m andV by V.
In summary, the offline POD-(M)DEIM is presented in Box 3. The procedure is as follows. First, N s sampling points in the parameter space Ω are generated, and for every point µ i , solve the corresponding FOM for u * ∈ [0, u max ], i.e., the entire loading path. Snapshots of the nodal displacement vector, nodal damage vector, external force, and global matrices K aa and K dd for each load increment (designated by t j ) are stored. These snapshots are next used to build the bases for the POD, i.e., V and V, and Φ f , Φ m , Φm for POD-(M)DEIM. Finally, the DEIM (Box 2) is applied to Φ f , Φ m , Φm to obtain the reduced matrices and vectors and the reduced mesh. Recall that m f , m a , mā are the number of basis vectors for the force vector, matrices K aa and K dd , respectively.
Once all the bases are obtained and stored, online simulations of the brittle fracture problem for any µ ∈ Ω can be performed using the procedure given in Box 4. We used a k r,n+1 to denote the reduced-order nodal displacement vector at load increment n + 1 of AM iteration k. This notation also applies to the reduced nodal damage vector. As can be seen that one has to modify the FOM code to have a POD-(M)DEIM ROM code. The meshfree Kriging model, presented in the next section, is a non-intrusive technique where one does not modify the FOM code.

Remark 3.
For anisotropic PFMs, the displacement sub-problem is nonlinear and most often solved with the Newton-Raphson (NR) method. For each NR iteration, one has to solve K aa T δa = f ext − f int (a) for the displacement corrections. Therefore, in the POD-(M)DEIM, one simply applies the DEIM to the internal force vector f int (a) and MDEIM to the tangent matrix K aa T . This constitutes a small modification to our proposed POD algorithm. Thus, our POD ROM is inherently multi-dimensional. See [47] for an application of POD to nonlinear finite elements.
Compress the local snapshots to the global ones:

Meshfree Kriging Method
This section presents the Kriging model that we utilized, for which more details can be found in [53,54,60]. Generally, in order to apply the Kriging prediction model, one follows the algorithm in Figure 1. Basically, it also consists of two phases: the offline phase where data are collected by running the FOM and the predictor is built. In the online phase, the predictor is used to get the outputs for any given input. Note that this phase does not use any phase-field finite element code, resulting in a non-intrusive model (see Remark 4).  Consider an output of interest y = {y(x; µ)|µ ∈ Ω} that varies within a parameter space Ω. Kriging models can be obtained by assuming that output y(x; µ) can be described as a linear combination of a regression model and a stochastic process as follows [60]: where m is a regression model known as the deterministic trend, which globally approximates the parameter space. The regression model is a linear combination of p chosen functions f j , where β j are regression parameters and f j are regressor functions. The stochastic process Z, which creates the localized deviation of the parameter space, is assumed to have zero mean and variance function: where σ 2 is variance and R() is the correlation model with parameters θ. We now consider a set of N s design points denoted by µ d = {µ 1 , . . . , µ N s } to which the corresponding output are y = {y(x; µ 1 ), . . . , y(x; µ N s )}. The regressor F with F ij = f j (x; µ i ) is given by a vector of p regressor functions: By defining the correlation matrix R(θ) as the matrix of the stochastic process correlation between the i th and j th design points, we have: The vector of correlations between an untried point, µ, and the design points is defined as: The predicted estimates,ŷ(x; µ), of the response y(x; µ) at an untried point µ, are given by: whereβ is estimated from the data, using least squares regression: Once the regression model and the covariance function of a stochastic process have been determined, the prediction of y can be done by using Kriging models. According to the Design and Analysis of Computer Experiments (DACE) [60], the correlations are restricted in the form: where d is the dimension of the parameter space. The correlation parameters can be determined by minimizing the log-likelihood for θ as: where |R(θ)| is the determinant of the correlation matrix corresponding to the design points. Assuming a Gaussian process, the optimal coefficientsβ of the correlation function solve: The corresponding maximum likelihood estimate of the variance, σ 2 , is: The model template in the DACE toolbox (Design and Analysis of Computer Experiments) has been discussed in [61].

Remark 4.
As the Kriging model is non-intrusive in the sense that one can just use a PFM code as a black box to obtain the training data, it applies to any phase-field models, namely 2D/3D and isotropic/anisotropic models.

A Posteriori Error Estimations
The relative L 2 -norm and Root Mean Squared Error (RMSE) used to evaluate our reduced-order models are written as: and: whereâ is the approximate ROM solutions and a represents the true solutions (precisely the FOM solutions). Note thatâ = Va r , so both a andâ have the same length.

Numerical Examples
In this section, the performances of the two proposed ROMs are evaluated for the 1D problem stated in Section 2. Even though one can, in general, consider the variation of the geometry (L), material parameters (E 0 , G c , b), and loadings (u * ), we present, for simplicity, results for µ = (E 0 , G c ) ∈ Ω ⊂ R 2 . That is, we consider a traction bar of a fixed length L = 25, b = L/100, unit cross-section (units are deliberately left out here, given that they can be consistently chosen in any system), and that is subjected to the maximum load of u max = 37.5. Material constants E and G c are varying parameters. Previous studies have shown that such a length scale is small enough [23] (see Remark 5). The bar is uniformly discretized with 1000 linear elements (element size h < b/2), resulting in N = 1001 grid points, which produce converged results; see Appendix B for a convergence analysis. The load is applied via prescribed displacement u * ∈ [0, u max ] at the right-most node. The entire loading process is divided into N T equal steps. The quantities of interests are (i) the displacement field, strain field, and damage field along the bar at u max and (ii) the stress-strain curve, which is obtained from the load-displacement curve where load is the reaction force at the right-most node and displacement is the applied displacement u * . These quantities computed using an FOM serve as the output of interest y = y(x; µ) used to build the Kriging predictor. All simulations were carried out using an in-house code written in MATLAB.
Concretely, we considered the parameter space E 0 × G c = [1,10] × [1,10]. The offline calculation can face the "curse of dimensionality" if we generate the samples using a full-factorial experiment. For example, let N E 0 = N G c = 50 be the sample numbers of each variable; the computations require N E 0 × N G c = 2500 runs for evaluating the solution for every possible combination of every possible design value. Therefore, in this paper, we have used the Latin Hypercube Sampling method (LHS) [62] to generate statistically-optimal sampling points. The typical behaviour of this traction softening bar is as follows; see Figures 5 and 6. When the load is small, i.e., smaller than 0.8u max , the bar is still in the elastic regime, albeit with a small damage due to the lack of an elastic domain of the chosen PFM. A further increment in loading leads to the initiation of a crack (a point in this 1D case) in the centre of the bar. Note that we have not introduced any imperfections in the bar to trigger damage localization; see, e.g., [55]. The fact that the point of damage localization is the centre of the bar is due to the perfect symmetry of the problem. The strain and damage are now localized in a small region centred at this cracking point.

Remark 5.
We note that for the selected phase-field model (Section 2), the solution depends on the length scale b (see Appendix B), and one should consider its variation. However, if one adopts our length scale-insensitive PFM presented in [7], the result is independent of b.

Remark 6.
In order to make sure the ROMs can capture strain localization (or crack initiation), we have applied a large displacement of u max = 37.5 to make sure that for all considered sampling points (i.e., all pairs of (E 0 , G c ) used in the offline phase), damage localization happens at least once.
We first present the ROM results in Section 4.1. That is, we generated sampling points, built various ROMs (with different accuracies), and selected the best ROM, by evaluating the ROMs with the FOM for a given random (E 0 , G c ), to be used in online computations. The construction of the Kriging's model was discussed in Section 4.2. Finally, Section 4.3 presents the comparison between the POD ROM and Kriging model by evaluating them against the FOM for some random values (E 0 , G c ) ∈ Ω.

ROM Results
To illustrate the performance of the POD-(M)DEIM approach, we generated N s = 60 sample points in Ω using the LHS, and built the ROM using the algorithm given in Box 3. This number of sampling points was based on our experiences with ROM; see for instance [42]. Precisely five ROMs, cf. Table 1, have been built to obtain the (M)DEIM bases of K aa , K dd , and f. Figure 2 shows the snapshot spectrum and the cut-off lines corresponding to the POD = 10 −9 of these matrices and vector. Concretely, we have used Step 1 in the algorithm in Box A1 (Appendix A) to get the snapshot spectrum and Step 2 to get the number of required basis vectors. Table 1     Now, for each (M)DEIM in ROM i , we built the POD basis of the solution's snapshot (displacement and damage fields) with the assumption that the number of POD bases for the displacement and and that for the damage field were equal, i.e., m =m. A number of cases corresponding to m ∈ {2, 100} were considered. The POD-(M)DEIM was then validated against the FOM and pure POD approach for a random pair of (E 0 , G c ). Figure 3 presents the L 2 -norm relative error (L 2 f for the reaction force) of ROM i , with i = 1, . . . , 5 and pure POD in comparison with FOM. In terms of accuracy, the pure POD performed excellently in the range of (20 − 40) POD basis. When more bases were used, the accuracy was reduced due to the fact that more noise and error were added into the model. Meanwhile, the accuracy of ROM i increased when the captured energy increased. However, the accuracy was not much different when POD was in the range of 10 −9 and 10 −10 . Furthermore, the performance of ROM i was much faster than the pure POD. We note that the performance (all of the computations were performed on a PC with Intel(R) Core(TM) i7-6820HQ CPU @ 2.70 GHz 2.7 GHz, RAM 8.0 GB (64-bit operating system, x64-based processor)), in terms of both accuracy and speed, of ROM 4 and ROM 5 was better than ROM 1 , ROM 2 , and ROM 3 . The reason here is that the AM solver (cf. Box 3)

Kriging Results
In order to verify the Kriging model against the sampling numbers, we generated several samples using LHS with N s from 10-5000, ran the FOM, and collected the output of interest, then built the Kriging model as discussed in Section 3.2; see also Figure 1. The second-order polynomial was chosen as a regression model, while the Gaussian function served as a correlation model. After that, a random pair of (E 0 , G c ) was selected to test the Kriging's predictor. Figure 4 shows the relative L 2 -norm and RMSE of the reaction force, the construction, and the prediction time, respectively. The relative error reduced when the number of samples increased (see Remark 7 for an elaboration on this error behaviour). We note here that the L 2 -norm error and the RMSE produced similar values. The prediction time was generally relatively constant; however, the model's construction time was increased when the sampling number increased, for example, t build = 0.1025 s for N s = 50 samples; however, it increased rapidly to t build = 1.17 × 10 3 s for N s = 4000. From here, to avoid over-fitting, we can use N s = 4000 (instead of N s = 1000) as our Kriging model to compare with ROM. Please note here that although t build was considered large for N s = 4000, it was computed one time at the offline phase, so it did not affect the performance of the prediction.

Remark 7.
Actually, the relative error in Figure 4 reduced significantly in the interval [1,1000] of samples, and then, it fluctuated, but did not go down to zero. From the data analysis viewpoint, the more data, the greater the accuracy of the model. However, it seemed to not be the case here. This was probably due to the discontinuity of the data (damage localization), and the Kriging model could not improve the accuracy around the discontinuity, although more data were introduced.

ROM vs. Kriging
A random sample test with N test = 10 was generated to compute the relative errors of the ROM and Kriging prediction in comparison with the FOM. The averaged errors and computational time of each model are recorded and given in Table 2. Although the ROM provided more accurate solutions than Kriging (for example, relative error L 2 f of 10 −4 vs. 10 −3 ), the ROM online phase was much slower than the Kriging model. Concretely, ROM's speedup factor, compared with FOM, was approximately eight-times, while the Kriging's was approximately 1100. Table 2 We now present the mechanical behaviour of this traction bar for the case (E 0 , G c ) = (3.705, 4.332). The response of the bar, obtained using the FOM, in terms of the load-displacement curve is given in Figure 5, where F is the reaction force at the right-most node. As this phase-field model lacked an elastic domain, the pre-peak behaviour was not linear, as damage was non-zero immediately upon load application. When the peak was reached, the bar was suddenly broken into two parts reflected by a sharp drop in the load-displacement curve. The evolution in time of the displacement, damage, and strain field is shown in Figure 6 for three time instances: 0.4u max , 0.8u max , and 0.9u max as marked in Figure 5. Evidently, there was a strong localization of damage and strain at the middle of the bar. It is interesting to see whether ROM solutions can capture this strain localization phenomenon.

Conclusions
Within the context of parametrized brittle fracture mechanics, we have presented two classes of reduced-order phase-field models. They can be used to carry out a very large number of computations required in different situations efficiently, ranging from parameter sensitivity analysis for brittle fracture problems, parameter estimation for phase-field brittle fracture models, to fracture constrained optimization problems. Our first ROM was a Proper Orthogonal Decomposition (POD)-based projection method that utilized the Discrete Empirical Interpolation Method (DEIM) and the Matrix Discrete Empirical Interpolation Method ((M)DEIM) to approximate the nonlinear vectors and matrices. The second was a meshfree Kriging model. For one-dimensional problems where Young's modulus and fracture energy vary, the POD-(M)DEIM ROM can speed up the online computations eight-times, whereas the Kriging model's speed up factor was 1100, albeit with a slightly lower accuracy. The greatest difference between the POD-(M)DEIM and Kriging model was the non-intrusive nature of the latter in the sense that one does not have to modify the full-order model code. Needless to say, all these computational savings were obtained with extensive offline computations using the full model.
Even though our reduced-order models were applied for a one-dimensional bar, we have shown that they are inherently multi-dimensional in nature. However, further works are required for 2D/3D fracture problems. It is possible that one might need to use more sampling points to capture different crack patterns. Additionally, the proposed models suffered from the following limitations

•
They cannot be used for extrapolation, i.e., when the parameters are out of the bounds of the considered parameter space; • The load has not been parametrized. That is, the maximum prescribed displacement u max is fixed. • The Kriging model resulted in oscillations around the damage localization point.
We note that the first limitation is inherent to any interpolation-based methods and might be tackled using machine learning methods. It is straightforward to overcome the second issue by just building a POD for u max . It is, however, very difficult to parametrize loadings in 2D and 3D. As far as the third issue is concerned, we anticipate that deep learning techniques might be helpful. For higher dimensional problems, it can become difficult. We are pursuing these paths and hope to publish them in the near future.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. Proper Orthogonal Decomposition Algorithm
For completeness, the POD algorithm to find the basis vectors directly from the snapshot matrix X ∈ R n×n s is given in Box A1. One can input either the desired number of basic vectors (m) or the tolerance POD . Note that n can be N for the displacement/damage/force snapshot or N 2 for the matrix snapshots (K aa , K dd ). n s is the number of snapshots. The total energy and capturing energy are denoted by a and b, respectively. c This is to get the first m eigenvalues of Σ.

Appendix B. Exact Solutions
For the model given in Section 2, one can derive the exact solution for the homogeneous stress state, which is given by [23,55]: where σ c denotes the maximum stress, which depends on b. Note that the homogeneous solution is valid up to the point of damage localization or the peak load in the load-displacement curve. We present mesh convergence analysis to justify the utilization of 1000 elements (or 1001 nodes) in all of our analyses. Three meshes consisting of 100, 500, and 1000 elements corresponding to the cases h = b, h = b/5, and h = b/10 were considered, where h denotes the element size. The results for the stress-strain are given in Figure A1a and the results for the damage profile in Figure A1b. As can be seen, 500 elements and 1000 elements yielded identical stress-strain data (thus, the same peak load), and 1000 elements were required to capture the cusp of the damage profile d(x) = exp(−|x|/b) [6]. Therefore, we denote the results obtained with 1000 elements' converged solutions and using this mesh to build the ROMs. Doing so, the error is due to ROM only, and thus, it is easier to understand the ROM behaviour.