Recursive multilevel trust region method with application to fully monolithic phase-field models of brittle fracture

TL;DR

This paper introduces a recursive multilevel trust region method to efficiently solve the non-convex minimization problems in phase-field models of brittle fracture, improving simulation of crack propagation in complex materials.

Contribution

It proposes a novel RMTR method with level-dependent objective functions tailored for phase-field fracture models, enhancing convergence and computational efficiency.

Findings

Method demonstrates improved convergence in 3D simulations.

Level-dependent objectives effectively combine fine and coarse mesh information.

Numerical examples validate the method's efficiency and robustness.

Abstract

The simulation of crack initiation and propagation in an elastic material is difficult, as crack paths with complex topologies have to be resolved. Phase-field approach allows to simulate crack behavior by circumventing the need to explicitly model crack paths. However, the underlying mathematical model gives rise to a non-convex constrained minimization problem. In this work, we propose a recursive multilevel trust region (RMTR) method to efficiently solve such a minimization problem. The RMTR method combines the global convergence property of the trust region method and the optimality of the multilevel method. The solution process is accelerated by employing level dependent objective functions, minimization of which provides correction to the original/fine-level problem. In the context of the phase-field fracture approach, it is challenging to design efficient level dependent objective functions as the underlying mathematical model relies on the mesh dependent parameters. We introduce level dependent objective functions that combine fine level description of the crack path with the coarse level discretization. The overall performance and the convergence properties of the proposed RMTR method are investigated by means of several numerical examples in three dimensions.

Tables

Table 1. Table 1: Setup of parameters used in numerical experiments.

Material parameter	Test case
	Fracture modes	Mixed-mode	Conchoidal fracture	Pressurized fracture
$\lambda \left(\text{GPa}\right)$	$12.1$	8.8	100	$12$
$\mu \left(\text{GPa}\right)$	$7.9$	1.3	100	$8$
$G_c\left(\text{GPa}\right)$	$5\times {10}^{-4}$	$1.1\times {10}^{-4}$	$1\times {10}^{-3}$	$1\times {10}^{-3}$

Table 2. Table 2: Choice of parameters used inside of TR/RMTR algorithms.

Parameter	${\eta }_1$	${\eta }_2$	${\gamma }_1$	${\gamma }_2$	${\Delta }_{0,0}^L$	${\mu }_1$	${\mu }_2$
Value	$0.1$	$0.75$	$0.5$	$2.0$	$1$	1	1

Table 3. Table 3: Speedup of RMTR method configured with the solution dependent 2nd order model with respect to alternative variants, namely 1st order consistency, 2nd order consistency and Galerkin model. Study performed on three examples: Mixed mode, Conchoidal fracture, and Pressurized fracture.

	Mixed mode	Conchoidal fracture	Pressurized fracture
1st order	$7.42$	$1.36$	$27.86$
2nd order	$2.13$	$1.12$	$2.12$
Galerkin	$2.65$	$1.30$	$5.33$

Table 4. Table 4: Number of dofs used for testing three fracture modes. The number of dofs is displayed with respect to the different refinement level. One should note, that the employed length-scale parameter l s subscript 𝑙 𝑠 l_{s} also changes accordingly.

Level	1	2	3	4	5	6
# dofs	$480$	$2,484$	$15,300$	$105,732$	$782,340$	$6,010,884$
$l_s$	$0.2$	$0.1$	$0.05$	$0.025$	$0.0125$	$0.00625$

Equations

$$\Psi (\mathbf{u}):=\underset{\text{elastic energy}}{\underbrace{{\int }_{\Omega }{\psi }_e(\mathbf{\epsilon }(\mathbf{u}))d\Omega }}+\underset{\text{fracture energy}}{\underbrace{{\int }_{\Gamma }\text{pazocal}{G}_{c}d\Gamma .}}$$

$${\int }_{\Gamma }\text{pazocal}{G}_{c}d\Gamma \approx \underset{\text{reaction term}}{\underbrace{{\int }_{\Omega }\frac{\text{pazocal}{G}_c}{2l_s}{c}^{2}d\Omega }}+\underset{\text{diffusion term}}{\underbrace{{\int }_{\Omega }\frac{\text{pazocal}{G}_c{l}_s}{2}|\nabla c{|}^{2}d\Omega .}}$$

$$\displaystyle\Psi({\mathbf{u}},c):=\underbrace{\int_{\Omega}[g(c)(1-k)+k]\psi_{e}^{+}(\bm{\varepsilon}({\mathbf{u}}))+\psi_{e}^{-}(\bm{\varepsilon}({\mathbf{u}}))\ d\Omega}_{\text{elastic energy}}+\underbrace{\int_{\Omega}\pazocal{G}_{c}\Big{[}\frac{1}{2l_{s}}c^{2}+\frac{l_{s}}{2}|\nabla c|^{2}\Big{]}\ d\Omega,}_{\text{fracture energy}}$$

$$(\mathbf{u},c)\in \text{argmin}\Psi (\mathbf{u},c)$$

$$\text{subject to}{\partial }_{t}c\ge 0,$$

$$\mathbf{U}\in \text{argmin}\Psi (\mathbf{U})$$

$$\text{subject to}\mathbf{U}\ge {\mathbf{U}}^{\mathbf{t}-1}\text{on}\Phi ,$$

$$\begin{array}[]{ll}\underset{{\mathbf{x}}\in{\mathbb{R}}^{n}}{\text{minimize}}&f({\mathbf{x}})\\[8.61108pt] \text{subject to}&{\mathbf{l}}\leq{\mathbf{x}}.\end{array}$$

$$m_i({\mathbf{p}}_i)=f_i+⟨{\mathbf{g}}_i,{\mathbf{p}}_i⟩+\frac{1}{2}⟨{\mathbf{p}}_i,{\mathbf{H}}_i{\mathbf{p}}_i⟩,$$

$$\text{pazocal}{B}_i:=\{{\mathbf{x}}_i+\mathbf{p}\in {\mathbb{R}}^n|\Vert \mathbf{p}{\Vert }_{\infty }\le {\Delta }_i\}.$$

$$\text{pazocal}{W}_i=\text{pazocal}{B}_i\cap \text{pazocal}F.$$

$$\begin{array}[]{ll}\underset{{\mathbf{p}}\in{\mathbb{R}}^{n}}{\text{minimize}}&m_{i}({\mathbf{p}})\\[8.61108pt] \text{such that }&{\mathbf{x}}_{i}+{\mathbf{p}}\in\pazocal{W}_{i}.\\[8.61108pt] \end{array}$$

$${\rho }_i=\frac{f({\mathbf{x}}_i)-f({\mathbf{x}}_i+{\mathbf{p}}_i)}{m_i(0)-m_i({\mathbf{p}}_i)}=\frac{\text{actual reduction}}{\text{predicted reduction}}.$$

$${\mathbf{x}}_{\nu ,0}^{l-1}={\mathbf{P}}_l^{l-1}{\mathbf{x}}_{\nu ,{\mu }_1}^l.$$

$$\begin{array}[]{ll}\underset{{\mathbf{p}}^{l-1}_{\nu}\in{\mathbb{R}}^{n^{l-1}}}{\text{minimize}}&\mathscr{h}_{\nu}^{l-1}({\mathbf{x}}_{\nu,0}^{l-1}+{\mathbf{p}}^{l-1}_{\nu})\\[8.61108pt] \text{such that}&{\mathbf{x}}_{\nu,0}^{l-1}+{\mathbf{p}}^{l-1}_{\nu}\in\pazocal{L}_{\nu}^{l}.\end{array}$$

$$\underset{{\mathbf{p}}_{\nu ,i}^l\in {\mathbb{R}}^{n^l}}{\mathrm{minimize}}$$

$${\mathbf{x}}_{\nu ,i}^l+{\mathbf{p}}_{\nu ,i}^l\in \text{pazocal}{W}_{\nu ,i}^l,$$

$${\mathbf{p}}_{\nu ,{\mu }_1+1}^l={\mathbf{I}}_{l-1}^l({\mathbf{x}}_{c,{\mu }^{l-1}}^{l-1}-{\mathbf{x}}_{\nu ,0}^{l-1}).$$

$$h_{\nu ,{\mu }_1+1}^l<h_{\nu ,{\mu }_1}^l.$$

$${\rho }_{\nu ,{\mu }_1+1}^l=\frac{h_{\nu }^l({\mathbf{x}}_{\nu ,{\mu }_1}^l)-h_{\nu }^l({\mathbf{x}}_{\nu ,{\mu }_1}^l+{\mathbf{p}}_{\nu ,{\mu }_1}^l)}{h_{\nu }^{l-1}({\mathbf{P}}_l^{l-1}{\mathbf{x}}_{\nu ,{\mu }_1}^l)-h_{\nu +1}^{l-1}({\mathbf{x}}_{\nu ,{\mu }^{l-1}}^{l-1})}=\frac{\text{fine level reduction}}{\text{coarse level reduction}}.$$

$$h_{\nu }^l({\mathbf{x}}_{\nu }^l):=\underset{\text{coarse model}}{\underbrace{f_{\nu }^l({\mathbf{x}}_{\nu }^l)}}+\underset{\text{1st order coupling}}{\underbrace{⟨\delta {\mathbf{g}}_{\nu }^l,{\mathbf{x}}_{\nu }^l-{\mathbf{x}}_{\nu ,0}^l⟩,}}$$

$$\delta {\mathbf{g}}_{\nu }^l:=\{\begin{array}{ll}{\mathbf{R}}_{l+1}^l{\mathbf{g}}_{\nu ,{\mu }_1}^{l+1}-\nabla f_{\nu }^l({\mathbf{x}}_{\nu ,0}^l)& \text{if}l<L\\ 0& \text{if}l=L.\end{array}$$

$$h_{\nu }^l({\mathbf{x}}_{\nu }^l):=⟨{\mathbf{R}}_{l+1}^l{\mathbf{g}}_{\nu ,{\mu }_1}^{l+1},{\mathbf{x}}_{\nu }^l-{\mathbf{x}}_{\nu ,0}^l⟩+\frac{1}{2}⟨{\mathbf{x}}_{\nu }^l-{\mathbf{x}}_{\nu ,0}^l,({\mathbf{R}}_{l+1}^l{\mathbf{H}}_{\nu ,{\mu }_1}^{l+1}{\mathbf{I}}_l^{l+1})({\mathbf{x}}_{\nu }^l-{\mathbf{x}}_{\nu ,0}^l)⟩,$$

$$h_{\nu }^l({\mathbf{x}}_{\nu }^l):=\underset{\text{coarse model}}{\underbrace{f_{\nu }^l({\mathbf{x}}_{\nu }^l)}}+\underset{\text{1st order coupling}}{\underbrace{⟨\delta {\mathbf{g}}_{\nu }^l,{\mathbf{x}}_{\nu }^l-{\mathbf{x}}_{\nu ,0}^l⟩}}+\underset{\text{2nd order coupling}}{\underbrace{\frac{1}{2}⟨{\mathbf{x}}_{\nu }^l-{\mathbf{x}}_{\nu ,0}^l,\delta \mathbf{H}({\mathbf{x}}_{\nu }^l-{\mathbf{x}}_{\nu ,0}^l)⟩.}}$$

$$\delta {\mathbf{H}}_{\nu }^l:=\{\begin{array}{ll}{\mathbf{R}}_{l+1}^l{\mathbf{H}}_{\nu ,{\mu }_1}^{l+1}{\mathbf{I}}_{l+1}^l-{\nabla }^2{f}_{\nu }^l({\mathbf{x}}_{\nu ,0}^l)& \text{if}l<L\\ 0& \text{if}l=L.\end{array}$$

$$h_{\nu }^l({\mathbf{x}}_{\nu }^l):=\underset{\text{modified coarse model}}{\underbrace{{\tilde{f}}_{\nu }^l({\mathbf{x}}_{\nu }^l)}}+\underset{\text{1st order coupling}}{\underbrace{⟨\delta {\mathbf{g}}_{\nu }^l,{\mathbf{x}}_{\nu }^l-{\mathbf{x}}_{\nu ,0}^l⟩}}+{\chi }_1(\mathbf{c},ϵ)\underset{\text{2nd order coupling}}{\underbrace{\frac{1}{2}⟨{\mathbf{x}}_{\nu }^l-{\mathbf{x}}_{\nu ,0}^l,\delta \mathbf{H}({\mathbf{x}}_{\nu }^l-{\mathbf{x}}_{\nu ,0}^l)⟩.}}$$

$${\chi }_1(\mathbf{c},ϵ):=\{\begin{array}{ll}0& \text{if}\underset{i=0,\cdots ,n^L}{\max (c_i)}>ϵ\\ 1& \text{otherwise},\end{array}$$

$$\displaystyle\tilde{\psi}({\mathbf{u}},c):=\int_{\Omega}\Big{(}[g(c)+k]\psi_{e}^{+}(\bm{\varepsilon}({\mathbf{u}}))+\psi_{e}^{-}(\bm{\varepsilon}({\mathbf{u}}))\Big{)}\ d\Omega+\chi_{1}\int_{\Omega}\Big{(}\pazocal{G}_{c}\Big{[}\frac{1}{2l_{s}}c^{2}+\frac{l_{s}}{2}|\nabla c|^{2}\Big{]}\Big{)}\ d\Omega,$$

$$(\Pi (m_h),q_h{)}_{L_2(\text{pazocal}I)}=(s_h,q_h{)}_{L_2(\text{pazocal}I)}\forall q_h\in \text{pazocal}{X}^q,$$

$${\int }_{\text{pazocal}I}(m_h-\Pi (m_h))q_{h}d\mathbf{x}={\int }_{\text{pazocal}I}(m_h-s_h)q_{h}d\mathbf{x}=0\forall q_h\in \text{pazocal}{X}^q.$$

$$\sum _{i=1}^{n^m}{m}_i{\int }_{\text{pazocal}I}{N}_i^m{N}_k^{q}d\mathbf{x}=\sum _{j=1}^{n^s}{s}_j{\int }_{\text{pazocal}I}{N}_j^s{N}_k^{q}d\mathbf{x}\text{for}k\in \{1,\dots ,n^q\},$$

$$\mathbf{B}\mathbf{m}=\mathbf{D}\mathbf{s},$$

Full text

Keywords: phase-field fracture, monolithic scheme, trust-region, multilevel, non-convex functional

1 Introduction

Predicting damage, crack and fragmentation patterns is a long-lasting challenge in computational mechanics. The phase-field approach to fracture addresses this task in an elegant way and has therefore become very popular. The basic idea behind the phase-field approach is to regularize sharp crack interfaces in the structure by means of smeared, volumetric crack zones. This removes the need to represent discontinuities and therefore avoids possible remeshing. The volumetric approximation of a crack surface is based on the Ambrosio-Torelli functional [1], which employs an auxiliary damage variable. This damage variable acts as an indicator, which marks the state of the material from intact to fully broken. The transition between the broken and intact part of the domain is modeled in diffused manner. The size of the transition zone is controlled by the mesh dependent length-scale parameter.

Several approaches to phase-field fracture can be found in the literature. In the physics community, evolution equations are derived by adapting the phase transition formalism of Landau and Ginzburg [2]. In the mechanics community, phase-field fracture approach is based on the variational formulation of the brittle fracture developed by Francfort and Marigo [3], which provides an extension to classical Griffiths theory of fracture [4]. The first numerical implementation of the variational phase-field approach is presented by Bourdin et al. in [5]. Miehe et al. [6, 7] modify the underlying mathematical model by introducing rate-independent formulation that ensures the local growth of the damage variable. Their formulation also introduces anisotropic energy functional, which allows for the degradation of the tensile energy with increased damage. Those modifications result into a thermodynamically consistent phase field formulation of the brittle fracture. Since then, the phase-field approach has become popular and was extended in various directions, such as dynamics [8] [9], large deformations [10] [11] [12], cohesive fracture [13] [14], shells and plates [15] [16], and used within several multiphysics applications [17] [18], [19] [20].

An important feature of any phase-field fracture simulator is to ensure crack irreversibility [6], which introduces a lower bound on the phase field variable. The resulting constraint can be treated directly in the solution strategy, for example by using active set [21], or semi-smooth Newton’s methods [22]. In order to simplify the solution process, alternatives, built on the modified energy, or weak formulations, are widely considered in the literature. For example, Bourdin [5] transforms inequality constraints into equality by using penalization term from [23]. A strategy based on augmented Lagrangian penalization is introduced by Wheeler et. al in [24]. Miehe et. al [7] handle the irreversibility by introducing so-called history variable, which replaces maximal tensile energy accumulated within the loading process. Disposing the inequality constraints usually causes deterioration of the conditioning of the Hessian matrix [25], or leads to inconsistencies between original energy functional and the modified weak form.

In this work, the phase-field formulation of fracture is based on the variational framework. The fracture problem is modeled by using an energy functional, whose minimizer corresponds to the associated solution. Algorithmically, the minimization process can be realized by using monolithic or staggered solution scheme. In the monolithic approach, the energy functional is minimized simultaneously for displacement and phase field [13]. This is numerically challenging, as the energy functional of the coupled problem is non-convex. In contrast, using a staggered solution scheme, the energy functional is minimized individually for the displacement and phase field, which gives rise to two convex minimization problems. This makes staggered solution scheme very robust and popular in the engineering community. However, Gerasimov et. al show in [26], that if convergent, the monolithic scheme can be faster than the staggered scheme since the staggered scheme tends to underestimate crack speed. Therefore, it is an important and challenging task to design efficient and reliable solution strategies for the fully monolithic scheme.

To our knowledge, there are only a few published results concerning the monolithic framework. These advances were accomplished either by modifying the underlying weak formulations or by applying some globalization strategy during minimization process. Vignollet et. al [27] augment the weak form by using a dissipation-based arc-length procedure [28]. Heister et. al apply "convexification" trick, where the phase field variable in the critical terms is replaced by its extrapolation in time [29]. Gerasimov et. al employ a line-search based Newton’s solver, which takes into account negative curvature [26]. A globalization of Newton’s method, based on error oriented approach from [30], is successfully applied by Wick in [31]. In [32], Wick employs so-called modified Newton’s method, which is based on dynamical switching between full Newton’s and modified Newton’s steps. The modified Newton’s steps are based on an altered Hessian matrix, where coupling terms are scaled in an appropriate manner.

All aforementioned globalizations of Newton’s method require the solution of a system of linear equations in order to determine step direction. Solving those linear systems can become computationally demanding as they are usually very large and ill-conditioned, see Section 2. An ideal choice of the linear solution strategy is a multilevel method due to its optimal complexity and scalability. Linear multilevel methods have been developed for more than 50 years and are used extensively in many areas. An introduction to multigrid methods can be found for example in [33, 34]. In the context of the phase field fracture models, an algebraic multilevel method is applied by Farrell and Maurini in [35]. The authors of [35] employ smoothed aggregation in order to precondition systems related to the displacement and the phase-field subproblems. Also within the staggered solution scheme, the geometric multigrid method is applied by Bilgen et. al in [36]. To the best of our knowledge, the multilevel solution strategy has not been applied yet to the coupled systems arising in the monolithic framework.

In this work, we present a globally convergent and fully-nonlinear multilevel solution method, for the phase field fracture problems originating in the monolithic framework. We employ a recursive multilevel trust region method (RMTR) [37, 38], which combines the trust region globalization strategy with the nonlinear multilevel framework. As other nonlinear multilevel methods, for instance FAS [39], or MG/OPT [40], the RMTR accelerates the solution process by employing level dependent objective functions. The minimization of these level dependent objective functions provides a coarse level corrections to the original/fine-level problem. Quality of the coarse level corrections is determined by the characteristics of the level dependent objective functions. A construction of an efficient level dependent objective functions for the phase-field fracture problems is a challenging task. This is due to the fact that the description of the fracture zones depends on the mesh-dependent length-scale parameter. The main contribution of this paper is a design of novel level dependent objective functions for the phase field fracture problems, see Section 4.

Employing recursive multilevel trust region method, with trust region radius defined by the $\ell_{\infty}$ norm, allows for a natural treatment of the variable bounds [41]. For this reason, the presented RMTR method is set-up in such a way, that all corrections satisfy the irreversibility condition. This is achieved by performing coarse level minimization with respect to a restricted irreversibility condition, see Section 3.

This paper is organized as follows: We start with a short review of phase-field fracture model, Section 2. In Section 3, we provide an introduction to the trust region and multilevel trust region methods. Section 5 discusses our choice of the transfer operators. Section 4 comments on how to build efficient level-dependent objective functions. The multilevel treatment of the irreversibility condition is addressed in Section 6. Numerical examples used to test the RMTR Algorithm will be introduced in Section 7. Finally, Section 8 demonstrates overall performance of the RMTR solution strategy. The presented work is summarized in Section 9.

2 Variational formulation of phase-field fracture

In this section, we briefly review the variational formulation of the second order phase-field model for brittle fracture. Our presentation focuses on quasi-static fracture problems and neglects inertia effects. A pseudo-time $t$ serves for indexing the deformation state of the loading process.

Let $\Omega\subset\mathbb{R}^{d},d=2,3$ be the computational domain with internal fracture surface $\Gamma\subset\mathbb{R}^{d-1}$. The computational domain $\Omega$ is assumed to have Lipschitz boundary $\partial\Omega$, which is decomposed into the Dirichlet $\partial\Omega_{D}$ and the Neumann $\partial\Omega_{N}=\partial\Omega\setminus\partial\Omega_{D}$ part. The variational approach to brittle fracture predicts the crack evolution by minimizing the total potential energy of the body $\Omega$, given by

[TABLE]

The energy functional (1) incorporates an elastic energy and a fracture energy. In this work, the elastic energy describes small deformations and utilizes the linearized strain tensor $\bm{\varepsilon}({\mathbf{u}})=0.5(\nabla{\mathbf{u}}+(\nabla{\mathbf{u}})^{T})$, where ${\mathbf{u}}$ is the displacement of a material point ${\mathbf{x}}\in\Omega$ at given pseudo-time $t$.

The term describing the fracture energy originates from the classical theory of brittle fracture [4, 42]. It states that the energy required to create a unit area of fracture surface equals the critical energy release $\pazocal{G}_{c}\in{\mathbb{R}}$. As a consequence, the fracture energy is defined as a contribution of $\pazocal{G}_{c}$ integrated over the fracture surface $\Gamma$. The presence of a fracture surface integral in (1) requires precise tracking of the crack surface paths and necessitates tedious remeshing, which leads into computationally exhausting algorithms.

2.1 Energy formulation of phase-field approach to fracture

The phase field approach to fracture [43] overcomes difficulties related to the remeshing of the crack surface $\Gamma$ by introducing the phase-field variable $c\in[0,1]$, which characterizes the material state from unbroken $c=0$ to fully broken $c=1$ and "damaged” in between. The regularization of the fracture is then performed by replacing fracture surface integral in (1) by its volumetric approximation

[TABLE]

Approximation (2) takes into account the phase-field $c$, its spatial gradient $\nabla c$ and regularization/length-scale parameter $l_{s}\in{\mathbb{R}}^{+}$. The reaction term in (2) corresponds to the area of the region, where material is damaged, e.g. $c\gg 0$. The diffusion term describes the area of the transition between the broken and intact part of the domain, where the gradient is high. The length-scale parameter $l_{s}$ controls the thickness of the damaged region by weighting the contributions of the reaction and the diffusion terms. Increasing $l_{s}$ makes the crack zone more diffused, while $l_{s}\rightarrow 0$ gives rise to a sharp crack surface $\Gamma$. The effect of the length scale parameter $l_{s}$ is demonstrated in Figure 2.

Taking into account (2), the total potential energy from (1) can be reformulated as follows

[TABLE]

where $g(c):=(1-c)^{2}$ describes the degradation of the stored energy due to evolving damage. It should be noted, that the degradation function $g(c)$ acts only on the tension part of the elastic energy. This is due to the fact, that in the reality, fracture develops only in the direction of tension [5]. The linear isotropic material models do not capture this feature and allow for fracture to occur in both tension and compression. In order to avoid unrealistic crack patterns, the anisotropic material models have to be employed, see for example [44], [45]. We employ an anisotropic model introduced by Miehe et. al in [6], where $\psi_{e}(\bm{\varepsilon}({\mathbf{u}}))$ is additively decomposed into $\psi_{e}^{+}(\bm{\varepsilon}({\mathbf{u}}))$ and $\psi_{e}^{-}(\bm{\varepsilon}({\mathbf{u}}))$ due to tension and compression, respectively. For more details about employed energy split, we refer interested reader to B.

Remark 1

For simplicity, formulation (3) neglects body forces and surface tractions. Artificial parameter $k\approx\mathscr{O}(\epsilon)$ in (3) circumvents a full degradation of energy at the completely broken state and helps to improve the condition number of the numerical system [1].

2.2 Minimization problem

At each loading step $t$, the displacement ${\mathbf{u}}$ and phase field $c$ are characterized as a minimizer of the following minimization problem:

Find $({\mathbf{u}},c)\in\ \bm{V}^{t}\times H^{1}(\Omega)$, such that

[TABLE]

where $\bm{V}^{t}:=\{{\mathbf{u}}\in\bm{H}^{1}(\Omega)\ |\ {\mathbf{u}}=\bm{g}^{t}\ \text{on}\ \ \partial\Omega_{D}\}$ is an admissible space for the displacement field. The symbol $H^{1}$ denotes the standard Sobolev space of weakly differentiable functions in $L^{2}$ with weak derivatives also in $L^{2}$. Since the minimization problem (4) is quasi-static, the time $t$ enters the above formulation only via the prescribed Dirichlet boundary condition $\bm{g}^{t}$ and the time-dependent irreversibility constraint $\partial_{t}c\geq 0$. This irreversibility constraint ensures positive evolution of the phase-field and prevents the crack from healing. Discretization of $\partial_{t}c\approx\frac{c^{t}-c^{t-1}}{\Delta t},$ with $\Delta t$ denoting time-step size allows for the reformulation of the minimization problem $\eqref{eq:min_problem_old}$ as follows:

Find $\bm{U}=({\mathbf{u}},c)\in\bm{V}^{t}\times H^{1}(\Omega)$ such that

[TABLE]

where $\bm{U^{t-1}}$ represents solution on the previous time-step and $\Phi=\bm{0}\times H^{1}(\Omega)$. In this way, the inequality constraint acts only the phase-field variable $c$.

2.2.1 Properties of the minimization problem

In order to solve the above minimization problem (5) numerically, we need to discretize underlying mathematical model. We perform spatial discretization of the problem by using finite-element method inside of the MOOSE framework [46]. We refer interested reader to A for the detailed formulation of the Euler-Lagrange equations.

Solving energy minimization problem (5) is numerically challenging for the following reasons:

• •

The functional (5) is non-convex and therefore its minimization admits many local minimizers. The most critical term is $[g(c)(1-k)+k]\psi_{e}^{+}(\bm{\varepsilon}({\mathbf{u}}))$ which is well-know to be non-convex in both variables ${\mathbf{u}}$ and $c$ [32]. The most popular solution strategy for solving nonlinear equations is Newton’s method, as it is simple, robust and quadratically convergent [25]. However, a convergence properties of the standard Newton’s method rely on a good initial guess. In the context of phase-field fracture, solver failures are typically reported in the case of brutal crack evolution and in the post-peak regime [26]. The convergence can be achieved by using small continuation steps [47] or by applying some globalization strategy, such as trust region [41], line-search methods [25], cubic regularizations [48] or affine invariant variants of the Newton method [30].

• •

Resolving the regularized crack surface $\Gamma_{l_{s}}$ demands a fine mesh with a high spatial resolution. In fact, the mesh size $h$ has to fulfill $h\leq\frac{l_{s}}{2}$. Here, the $l_{s}$ is the length-scale parameter from (2) controlling the width of the diffusive crack zone. As a consequence, we end up with large-scale coupled nonlinear systems of equations, which are computationally demanding to solve.

• •

Arising systems of equations might also be ill-conditioned as the presence of the crack leads to the strongly varying elastic stiffness. A suitable preconditioner, e.g. mutilevel [49], or domain decomposition methods [50] can to be employed to address this issue.

• •

The irreversibility condition defines a lower bound on the phase-field variable $c$. Therefore, any applied solution strategy, for example line-search method, has to take into account point-wise constraint. The popular choice include for instance active set methods [21], [22], interior point methods [51], or sequential quadratic programming methods [41].

This work addresses all aforementioned points by employing the recursive multilevel trust region algorithm (RMTR) [52, 38]. On one side, the multilevel framework addresses issues of the ill-conditioning and the large number of dofs. On other side, the trust region framework tackles the non-convexity of the energy functional. In addition, the trust region method with the trust region radius defined by the $\ell_{\infty}$ norm allows for a natural treatment of the point-wise constraints [41], see the Section 3.

3 Recursive Multilevel Trust Region

In this work, we propose to solve the energy minimization problem (5) with the recursive multilevel trust region method (RMTR) originally introduced in [37], [38]. Extended details about the method and its convergence properties can be also found in [52] and [53]. The variant of the RMTR method considered here is associated with the following optimization problem

[TABLE]

The objective function $f:{\mathbb{R}}^{n}\rightarrow{\mathbb{R}}$ in (6) denotes the energy functional (3). The solution vector ${\mathbf{x}}\in{\mathbb{R}}^{n}$ contains coefficients of the displacement and phase-field. The lower bound ${\mathbf{l}}\in{\mathbb{R}}^{n}$ originates from the crack irreversibility condition. One should note, that the constraint in (6) is defined component-wise, thus each component $({\mathbf{x}})_{k}$ of ${\mathbf{x}}$, where $k\in\{1,\dots,n\}$, has to satisfy $({\mathbf{l}})_{k}\leq({\mathbf{x}})_{k}$.

3.1 Trust region method

A trust region method (TR) is a step-size control strategy for minimizing nonlinear, convex/non-convex objective functions [41]. At each iteration $i$, the TR method approximates the objective function $f$ by a model $m_{i}$. The model $m_{i}$ is often constructed as a second order Taylor approximation of $f_{i}$ around the current iterate ${\mathbf{x}}_{i}\in{\mathbb{R}}^{n}$

[TABLE]

where we use the shortened notation $f_{i}=f({\mathbf{x}}_{i})$, ${\mathbf{g}}_{i}=\nabla f({\mathbf{x}}_{i})$ and ${\mathbf{H}}_{i}=\nabla^{2}f({\mathbf{x}}_{i})$. In this case, the update, i.e. the correction ${\mathbf{p}}_{i}$, is then obtained by minimizing the model $m_{i}$ within a certain region. This region is called the trust region and it is defined around current iterate ${\mathbf{x}_{i}}$. The size of the region is defined by the means of radius $\Delta_{i}>0$, which provides a constraint on the step size, e.g. $\|{\mathbf{p}}_{i}\|_{p}\leq\Delta_{i}$.

The choice of norm $\|\cdot\|_{p}$ defines shape of the trust region. In this work, we employ $\ell_{\infty}$ norm, e.g. the trust region is defined as

[TABLE]

As we will see later, defining trust region $\pazocal{B}_{i}$ as in (8) simplifies the treatment of point-wise constraint, e.g. irreversibility constraint. We can view (8) as point-wise constraints on step ${\mathbf{p}_{i}}$, thus $-\Delta_{i}\leq({\mathbf{p}_{i}})_{k}\leq\Delta_{i},\ \ \forall k\in\{1,\dots,n\}$.

In order to keep our iterates feasible throughout the computation, the minimization of model $m_{i}$ has to take place in the intersection of the feasible set $\pazocal{F}:=\{{\mathbf{x}}\in{\mathbb{R}}^{n}\ |\ {\mathbf{l}}\leq{\mathbf{x}}\}$ with the trust region $\pazocal{B}_{i}$. For this reason, we define the active working set $\pazocal{W}_{i}$ as follows

[TABLE]

Please note, that the working set $\pazocal{W}_{i}$ needs to be recomputed in each iteration, as the trust region $\pazocal{B}_{i}$ is iteration dependent. The evaluation of $\pazocal{W}_{i}$ is straightforward and can be performed component-wise.

At each iteration $i$, we obtain the step direction ${\mathbf{p}}_{i}$ by approximately solving the following constrained quadratic minimization problem

[TABLE]

The iterate obtained in this way (also called the trial point), ${\mathbf{x}}_{i}+{\mathbf{p}}_{i}$, is not taken immediately. Instead, the TR algorithm decides whether the trial point should be accepted and how to adjust the trust-region radius $\Delta_{i}$. To this aim, a trust-region ratio $\rho_{i}$, is defined as the quotient of actual reduction in $f_{i}$ and the reduction predicted by the model $m_{i}$. Given a step ${\mathbf{p}}_{i}$, $\rho_{i}$ is defined as

[TABLE]

The value of the trust region ratio $\rho_{i}$ close to 1 indicates good agreement between the model $m_{i}$ and the function $f$. In this case, it is safe to accept the trial step and expand the trust region radius. On the other hand, if value of $\rho_{i}$ is negative or close to 0, algorithm must reject trial step and shrink the trust region radius. Algorithm 1 summarizes the described process [25].

3.2 Nonlinear multilevel framework

Multilevel solution strategies take into account a hierarchy $\{\pazocal{X}^{l}\}_{l=0}^{L}$ of usually nested finite element, or approximation spaces. Here, we use the subscript $l$ to denote a given level and assume $L\geq 1$ to be the finest level. In the case of geometric multilevel methods, the hierarchy $\{\pazocal{X}^{l}\}_{l=0}^{L}$, with $\text{dim}(\pazocal{X}^{l})>\text{dim}(\pazocal{X}^{l-1})$ , is constructed from a hierarchy of meshes $\{\pazocal{T}^{l}\}_{l=0}^{L}$. For purpose of this work, the hierarchy of meshes $\{\pazocal{T}^{l}\}_{l=0}^{L}$ was obtained by uniformly refining an initial coarse mesh $\pazocal{T}^{0}$.

Transfer of the data between subsequent levels of multilevel hierarchy is ensured by transfer operators. In the presented work, we use three types of the transfer operators:

• •

The prolongation operator ${\mathbf{I}}_{l}^{l+1}:{\pazocal{X}}^{{l}}\rightarrow{\pazocal{X}}^{{l+1}}$ is designed to transfer primal variables, such as corrections from a coarse level to a finer level.

• •

The restriction operator ${\mathbf{R}}^{l}_{l+1}:\pazocal{X}^{{l+1}}\rightarrow\pazocal{X}^{{l}}$ transfers dual variables, e.g. gradients, from a fine level to a coarser level. In this work, the operators ${\mathbf{R}}^{l}_{l+1}$ and ${\mathbf{I}}_{l}^{l+1}$ are related by relation ${\mathbf{R}}^{l}_{l+1}=({\mathbf{I}}_{l}^{l+1})^{T}$.

• •

The projection operator ${\mathbf{P}}^{l}_{l+1}:\pazocal{X}^{{l+1}}\rightarrow\pazocal{X}^{{l}}$ transfers primal variables, for example the current iterate, from a fine level to a coarse level.

In our approach, all three transfer operators are computed based on pseudo-L2-projections, see Section 5.

The multilevel solution process is built as an interplay between smoothing and coarse grid correction. The smoothing step reduces high-frequency error related to each level. In the context of the trust region framework considered here, smoothing is associated with performing a small number of trust region iterations, where the arising quadratic models (7) are solved approximately. On the coarsest level, the minimization of coarse level objective function is again achieved by the trust region method, but constrained quadratic subproblems are solved accurately.

Minimization on each level $l$ is associated with level-dependent objective function $\mathscr{h}^{l}$. The definition of the level-dependent objective functions $\{\mathscr{h}^{l}\}_{l=0}^{L-1}$ is often based on the coarse level representation/discretization of the original model, denoted by $\{f^{l}\}_{l=0}^{L}$. A detailed discussion on how to construct the level-dependent objective functions $\{\mathscr{h}^{l}\}_{l=0}^{L-1}$ will be provided in Section 4.

3.3 RMTR algorithm

Our description of the RMTR, Algorithm 3, has the form of V-cycle and closely follows [38]. In the following, we use a triplet (cycle, level, iterate) to denote all computational quantities. For example the solution computed as the $i$-th iterate during the $\nu$-th cycle on level $l$ is denoted by ${\mathbf{x}}^{l}_{\nu,i}$.

Each V-cycle consists of nonlinear smoothing steps and a coarse level solve. Starting from the finest level, $l=L$, the RMTR algorithm performs $\mu_{1}\in\mathbb{N}$ pre-smoothing steps. This is done by a trust region solver, see Algorithm 1. The result of pre-smoothing is an intermediate iterate ${\mathbf{x}}^{l}_{\nu,\mu_{1}}$. The projection of this iterate to the next coarser level yields the initial guess on level $l-1$, i.e.

[TABLE]

Once we have reached the coarser level, we can perform the minimization of the level dependent objective function $h^{l-1}_{\nu}$. The result of this minimization, a coarse level correction, ${\mathbf{p}}_{\nu}^{l-1}={\mathbf{x}}^{l-1}_{\nu,\mu_{1}}-{\mathbf{x}}^{l-1}_{\nu,0}$, is afterwards prolongated back to the finer level. As a consequence, we have to guarantee that the the prolongated coarse level correction satisfies the following requirements:

It does not exceed the current fine level trust region radius, e.g. $\|{\mathbf{I}}_{l-1}^{l}{\mathbf{p}}^{l-1}_{\nu}\|_{\infty}\leq\Delta^{l}_{\nu,\mu_{1}}.$ 2. 2.

It does not violate the fine level irreversibility condition, e.g. ${\mathbf{I}}_{l-1}^{l}{\mathbf{p}}^{l-1}_{\nu}\in\pazocal{F}^{l}_{\nu}.$

We can meet both requirements by minimizing the coarse level objective function within a feasible set $\pazocal{L}_{\nu}^{l}$ as follows

[TABLE]

We will discuss in Section 6 how to construct set $\pazocal{L}_{\nu}^{l}$, such that $\|{\mathbf{I}}_{l-1}^{l}{\mathbf{p}}^{l-1}_{\nu}\|_{\infty}\leq\Delta^{l}_{\nu,\mu_{1}}$ and ${\mathbf{I}}_{l-1}^{l}{\mathbf{p}}^{l-1}_{\nu}\in\pazocal{F}^{l}_{\nu}$.

The minimization of $\mathscr{h}^{l}_{\nu}$ on a given level $l$ is performed by means of the trust region, Algorithm 1. The correction ${\mathbf{p}}_{\nu,i}^{l}$ is obtained by minimizing the quadratic model $m$ as

[TABLE]

where $\pazocal{W}_{\nu,i}^{l}$ is the working active set. As for the single level case, the set $\pazocal{W}_{\nu,i}^{l}$ is result of the intersection of the feasibility set $\pazocal{L}_{\nu}^{l}$ with the current trust region $\pazocal{B}_{c,i}^{l}:=\{{\mathbf{x}}_{\nu,i}^{l}+{\mathbf{p}}_{\nu,i}^{l}\in{\mathbb{R}}^{n^{l}}\ |\ \|{\mathbf{p}}_{\nu,i}^{l}\|_{\infty}\leq\Delta_{\nu,i}^{l}\}$, e.g. $\pazocal{W}_{\nu,i}^{l}:=\pazocal{L}_{\nu}^{l}\cap\pazocal{B}_{c,i}^{l}$.

3.3.1 Acceptance of the coarse level correction

After the approximate minimization on level $(l-1)$ is performed, we prolongate the resulting coarse level correction back to the finer level, i.e.

[TABLE]

The symbol $\mu^{l-1}$ in (16) denotes the number of all computed corrections on level $l-1$ during cycle $\nu$. The coarse level correction is not accepted immediately, but only if it provides a decrease in the fine level objective function; thus only if

[TABLE]

The quality of the prolongated coarse level correction ${\mathbf{p}}^{l}_{\nu,\mu_{1}+1}$ is determined by means of the multilevel trust region ratio $\rho^{l}_{\nu,\mu_{1}}$, which measures the agreement between the fine level and the coarse level reduction

[TABLE]

The correction ${\mathbf{p}}_{\nu,\mu_{1}+1}^{l}$ will be accepted only if $\rho_{\nu,\mu_{1}+1}^{l}>\eta_{1}$. This is in contrast to the classical nonlinear multilevel schemes, e.g. FAS [39], or NMG [34], where the coarse level correction is used independently of its quality. A similar strategy has been used in [54] for controlling global convergence.

Once the correction ${\mathbf{p}}_{\nu,\mu_{1}+1}^{l}$ has been accepted or rejected, the RMTR algorithm adjusts the trust region radius $\Delta_{\nu}^{l}$ accordingly, see Algorithm 2. The V-cycle completes with $\mu_{2}\in\mathbb{N}$ post-smoothing steps. The outlined process is performed recursively on each level, except on the coarsest one, where no recursion is invoked. On the finest level, the last iterate of cycle $\nu$ is taken over to the next cycle, $(\nu+1)$, as the initial iterate.

Termination of the RMTR algorithm is specified in terms of the criticality measure $\pazocal{E}$. Here, we follow [41] and define $\pazocal{E}({\mathbf{x}},f):=\|\mathscr{P}({\mathbf{x}}-\nabla f({\mathbf{x}}))-{\mathbf{x}}\|,$ where $\mathscr{P}$ is an orthogonal projection onto the feasible set. The details on how to compute $\mathscr{P}$ can be found in C.

4 Level dependent objective functions

On each level $l$, the RMTR Algorithm 3 minimizes some nonlinear objective function $\mathscr{h}_{\nu}^{l}$. The result of this minimization, the coarse level correction, is expected to provide decrease in the objective function on the subsequent finer level. The overall efficiency of the multilevel method relies on the quality of this prolongated coarse level correction. Our aim is to find coarse level model (level dependent objective function), which yields adequate coarse level correction. In this section, we review several models, which are often used in literature for constructing $\{\mathscr{h}_{\nu}^{l}\}_{l=0}^{L}$. Our discussion mostly focuses on the convergence properties and on the computational complexity of the given models. Later in this section, we propose an alternative model, which is specially designed for the phase-field fracture models considered in this work. All models discussed in this section assume $\mathscr{h}^{L}:=f^{L}$.

4.1 First order consistency

Classical nonlinear multilevel schemes, such as FAS [39], NMG [34], NMM [60] or MG/OPT [40], often employ a first order consistency model, where the level dependent function $\mathscr{h}^{l}:{\mathbb{R}}^{n^{l}}\rightarrow{\mathbb{R}}$ is defined as

[TABLE]

where

[TABLE]

Here, the first term, $f_{\nu}^{l}({\mathbf{x}}_{\nu}^{l})$, is the original model problem discretized on level $l$, thus $f_{\nu}^{l}$ is the energy functional $\Psi({\mathbf{u}},c)$ evaluated at the level $l$. The second term, $\langle\delta{\mathbf{g}}^{l}_{\nu},{\mathbf{x}}^{l}_{v}-{\mathbf{x}}^{l}_{\nu,0}\rangle$, couples the coarse level with the fine level. This is achieved by using the quantity $\delta{\mathbf{g}}^{l}_{\nu}$, which expresses the difference between the restricted fine level gradient and the initial coarse level gradient (evaluated at ${\mathbf{x}}^{l}_{\nu,0}={\mathbf{P}}^{l}_{l+1}{\mathbf{x}}^{l+1}_{\nu,\mu_{1}}$). The presence of this term enforces the relation ${\mathbf{g}}^{l-1}_{\nu,0}=\nabla\mathscr{h}_{\nu,0}^{l-1}({\mathbf{x}}^{l-1}_{\nu,0})={\mathbf{R}}^{l-1}_{l}{\mathbf{g}}_{\nu,\mu_{1}}^{l}$. We refer to the model (18) as first order consistency since the first step of the coarse level minimization process will be done in the direction of the negative restricted fine level gradient.

Model (18) is known to yield fast convergence, even with a poor initial guess, under the assumption that $h_{\nu}^{l}$ provides a good approximation of the error after pre-smoothing. If the coarse level model $h_{\nu}^{l}$ has poor approximation properties, then the multilevel method might produce an inadequate coarse level correction. This results in a slow convergence rate, even if the initial error is very low and nonlinearity is weak [61].

In terms of the phase-field fracture problems considered here, the model $h^{l}_{\nu}$ based on the $f^{l}_{\nu}$, the coarse level discretization of the energy functional $\Psi({\mathbf{u}},c)$, cannot be expected to provide good coarse level correction. This is due to the fact that volumetric approximation to the fracture surface integral (2) is determined by the mesh dependent length scale parameter $l_{s}$. Therefore, we are not able to represent the same fracture zones on the coarse level, as we can on the finest level. Figure 4 demonstrates this effect on two dimensional tension test. In particular, we consider a rectangular computational domain discretized with quadrilateral finite elements. The body is subjected to tensile loading, which results in a straight crack path inside the specimen. We perform the same experiment repeatedly while increasing the resolution of the mesh. In the same time, we also modify the length-scale parameter $l_{s}$ as follows $l_{s}=2h$, where $h$ represents mesh size. As we can see, the obtained crack paths differ for different resolution levels. Moreover, even larger inconsistencies between different refinement levels can be expected for more complex loading scenarios.

4.2 Galerkin approximation

Exploiting ideas from linear multilevel methods, the coarse level model $\mathscr{h}^{l}:{\mathbb{R}}^{n^{l}}\rightarrow{\mathbb{R}}$ can be created as a restricted version of the fine level model

[TABLE]

where ${\mathbf{R}}^{l}_{l+1}{\mathbf{g}}_{\nu,\mu_{1}}^{l+1}$ represents the restricted fine level gradient and ${\mathbf{R}}^{l}_{l+1}{\mathbf{H}}_{\nu,\mu_{1}}^{l+1}{\mathbf{I}}_{l}^{l+1}$ is the restricted fine level Hessian. In the unconstrained convex case, the formulation (20) is equivalent to applying a Galerkin linear multigrid on the associated Newton’s equation.

The overall performance of the method is determined by the accuracy of the fine level linearization, which depends on the properties of the Hessian and on the given initial guess. If the initial guess is in the local neighborhood of the solution, the method can be expected to achieve the quadratic convergence rate. However, for generic non-linear, non-convex functions, the 1st order consistency formulation (18), which tackles nonlinearity directly on the coarse level, yields usually better convergence rates [61]. From the computational point of view, the model (20) is simple to implement, as the multilevel solver does not require knowledge about coarse level discretization. Concerning the phase-field fracture simulations, the coarse level model (20) is able to represent fracture zones well, as it is built from the fine level data.

4.3 Second order consistency

Yavneh and Dardyk in [61] explore an idea of combining the first order consistency, e.g. (18), with the Galerkin approximation approach (20) by defining $\mathscr{h}^{l}:{\mathbb{R}}^{n^{l}}\rightarrow{\mathbb{R}}$ as follows

[TABLE]

Here, the term $\delta{\mathbf{g}}^{l}_{\nu}$ is as defined in (19). The quantity $\delta{\mathbf{H}}^{l}_{\nu}$ expresses the difference between restricted fine level Hessian and the initial coarse level Hessian

[TABLE]

The presence of this term enforces the first coarse level correction to be the restricted Newton’s direction. Subsequent coarse level corrections involve knowledge from the coarse level discretization of the energy functional. The multilevel method with $\{\mathscr{h}^{l}\}_{l=0}^{L}$ defined as in (21) is demonstrated to guarantee at least as good and usually better convergence rate than multilevel method with first order consistency or with the Galerkin approximation model, for the details see [61]. Despite its efficiency, the method is computationally more costly than previously discussed methods. This is due to the fact, that the model (21) requires knowledge about the coarse level discretization as well as the information about the matrix $\delta{\mathbf{H}}^{l}_{\nu}$. In addition, the product consisting of $\delta{\mathbf{H}}^{l}_{\nu}$ is used to update the coarse level gradient and the Hessian at each successful coarse level iteration.

4.4 Solution dependent second order consistency

In this section, we introduce a novel model in order to define of coarse level objective function $\mathscr{h}^{l}:{\mathbb{R}}^{n^{l}}\rightarrow{\mathbb{R}}$. The proposed model is basically variation of the (21) defined as following

[TABLE]

The two main ingredients of the formulation (23) are

• •

the indicator function $\chi_{1}$ is used to turn on and off the second order coupling term.

• •

the objective function ${f}_{\nu}^{l}({\mathbf{x}}_{\nu}^{l})$ is replaced by its modified variant $\tilde{f}_{\nu}^{l}({\mathbf{x}}_{\nu}^{l})$.

The first modification is justified by the fact that during the pre-stage phase of the simulation, the crack has not yet propagated. Thus, the coarse level discretization is a valid representation of the fine level model, with respect to the discretization error. In this case, the use of the first order consistency approach is sufficient and computationally cheaper. To this aim, we define indicator function

[TABLE]

which allows for turning on and off the second order consistency term. The symbol $c_{i}$ in (24) denotes the nodal coefficients of the fine level phase-field. A threshold $\epsilon\in{\mathbb{R}}$, $0\leq\epsilon<1$ defines a limit, after which we believe that the fine level model can not be represented well by the coarse level discretization. In our experience, $\epsilon\approx 0.85$ provides reasonable choice.

Our second modification consists of replacing the objective function ${f}_{\nu}^{l}({\mathbf{x}}_{\nu}^{l})$ by its modified variant $\tilde{f}_{\nu}^{l}({\mathbf{x}}_{\nu}^{l})$. As the staggered solution scheme [7], we split the total potential energy (3) into the elastic and the fracture energy. For the elastic part of the energy, we employ a coarse level discretization, since it can represent the fine level elastic energy well. On the other hand, for the fracture energy our coarse model uses only information restricted from the fine level, since the regularized fracture energy differs at the different discretization levels. To incorporate both of the above requirements, the modified objective function $\tilde{f}_{\nu}^{l}$ is defined as the modified energy functional

[TABLE]

where the indicator function $\chi_{1}$ is defined as in (24).

As it will be demonstrated in Section 8.1, the coarse level model (23), provides better coarse level error approximation when compared to the coarse level models (20), (18), (21). Conceptually similar approach was applied in the context of multilevel method for the frictional contact in [54], see also [60].

5 Transfer operators

The pseudo-$L^{2}$-projection belongs to the more generic group of mortar/variational transfer operators, which are commonly used in the context of non-conforming domain decomposition methods and contact problems [55]. The mortar projections allow for transferring the information of functional quantities from one approximation space to the another. Following the naming convention in the literature on the variational transfer, we transfer the data from the master space to the slave space. Therefore, in the following discussion, we use subscripts $m$ and $s$ to denote quantities related to the master and slave side, respectively.

5.1 Variational transfer

Let us consider meshes $\pazocal{T}^{m}$ and $\pazocal{T}^{s}$ with corresponding finite element spaces $\pazocal{X}^{m}$ and $\pazocal{X}^{s}$. The mortar projection $\Pi:\pazocal{X}^{m}\rightarrow\pazocal{X}^{s}$ maps a function from $\pazocal{X}^{m}$ to $\pazocal{X}^{s}$. For a function $m_{h}\in\pazocal{X}^{m}$, we are looking for $s_{h}=\Pi(m_{h})\in\pazocal{X}^{s}$, such that the following equality condition holds

[TABLE]

where $\pazocal{I}=\pazocal{T}^{m}\cap\pazocal{T}^{s}$. The space $\pazocal{X}^{q}$ in (26) denotes the multiplier space. Reformulating (26) results in the following weak equality condition

[TABLE]

Expressing $m^{h}\in\pazocal{X}^{m}$, $s^{h}\in\pazocal{X}^{s}$ and $q^{h}\in\pazocal{X}^{q}$ in terms of their respective basis, $\{N^{m}_{i}\}_{i=1}^{n^{m}}$, $\{N^{s}_{j}\}_{j=1}^{n^{s}}$, $\{N^{q}_{k}\}_{k=1}^{n^{q}}$, we obtain $m^{h}=\sum_{i=1}^{n^{m}}N^{m}_{i}m_{i}$, $s^{h}=\sum_{j=1}^{n^{s}}N^{s}_{j}s_{j}$, $q^{h}=\sum_{k=1}^{n^{q}}N^{q}_{k}q_{k}$, where $\{m_{i}\}_{i=1}^{n^{m}}$, $\{s_{j}\}_{j=1}^{n^{s}}$, $\{m_{k}\}_{k=1}^{n^{q}}$ are the coefficients. This allows us to reformulate equation (27) as

[TABLE]

which can be expressed by using vector-matrix notation as

[TABLE]

where the vectors ${\mathbf{m}}$ and $\mathbf{s}$ contain the coefficients $m_{i}$ and $s_{j}$. The components of the coupling matrix ${\mathbf{B}}$ are defined as $B_{k,i}=\int_{\pazocal{I}}N^{m}_{i}N^{q}_{k}\ d{\mathbf{x}}$, while the matrix ${\mathbf{D}}$ contains the entries $D_{k,j}=\int_{\pazocal{I}}N^{s}_{j}N^{q}_{k}\ d{\mathbf{x}}$. As a result, the discrete projection operator has the following algebraic representation

[TABLE]

5.2 Transfer operators via pseudo-$L^{2}$-projection

Different choices of the multiplier space $\pazocal{X}^{q}$ result in different transfer operators. For example, if we set $\pazocal{X}^{q}=\pazocal{X}^{s}$, the operator ${\mathbf{\Pi}}$ is nothing other than the well-known $L^{2}$ projection. In this particular case, the matrix ${\mathbf{D}}$ is the mass matrix on the slave space. For the purpose of this work, we consider the computationally cheaper alternative, e.g. pseudo-$L^{2}$ projection, where the multiplier space $\pazocal{X}^{q}$ is associated with the dual basis [56, 57]. The dual basis functions satisfy the following bi-orthogonality condition

[TABLE]

which gives rise to a diagonal matrix ${\mathbf{D}}$ in our setting. This is very convenient, as it makes an evaluation of the inverse in (30) computationally cheap. In addition, it produces a transfer operator with a sparse structure, which is not the case for the $L^{2}$-projection. The parallel implementation of the mortar projections is not a trivial task, however, a library MOONoLith [58] provides an efficient and scalable implementation.

As described in the Section 3.2, the RMTR Algorithm 3 makes use of three different transfer operators: prolongation, restriction, and projection. The prolongation operator ${\mathbf{I}}_{l}^{l+1}:{\pazocal{X}}^{{l}}\rightarrow{\pazocal{X}}^{{l+1}}$ transfers primal variables from a coarse level to the subsequent fine level. The assembly of the prolongation operator ${\mathbf{I}}_{l}^{l+1}$ is performed by the formula (30), with $\pazocal{X}^{m}=\pazocal{X}^{{l}}$ and $\pazocal{X}^{s}={\pazocal{X}}^{{l+1}}$. On the other hand, in order to assemble the projection operator ${\mathbf{P}}^{l}_{l+1}:\pazocal{X}^{{l+1}}\rightarrow\pazocal{X}^{{l}}$, we associate master and slave spaces as follows: $\pazocal{X}^{m}=\pazocal{X}^{{l+1}}$ and $\pazocal{X}^{s}={\pazocal{X}}^{{l}}$. The resulting operator ${\mathbf{P}}^{l}_{l+1}$ is then used to transfer iterates from the fine level to the coarse level. The restriction operator ${\mathbf{R}}^{l}_{l+1}:\pazocal{X}^{{l+1}}\rightarrow\pazocal{X}^{{l}}$, which is designed for transferring the dual variables from the fine level to the coarse is not assembled by using (30). Instead, it is obtained as the adjoint of the prolongation operator, thus ${\mathbf{R}}^{l}_{l+1}=({\mathbf{I}}_{l}^{l+1})^{T}$.

It is important to realize, that ${\mathbf{P}}^{l}_{l+1}\neq({\mathbf{I}}_{l}^{l+1})^{T}$. As noted by Groß and Krause in [38], the restriction operator ${\mathbf{R}}^{l}_{l+1}=({\mathbf{I}}_{l}^{l+1})^{T}$ is not suitable for transferring the primal variables from the fine level to the coarse level. This is due to the fact, that the restriction operator ${\mathbf{R}}^{l}_{l+1}$ is designed as a dual operator. The numerical results reported in [38] also demonstrate that the usage of the projection operator ${\mathbf{P}}^{l}_{l+1}$ based on $L^{2}$-projection accelerates the convergence behavior of the RMTR method. The detailed discussion about how to construct suitable operator ${\mathbf{P}}^{l}_{l+1}$ can be also found in [59].

Remark 2

We note, that we refer to the operator ${\mathbf{P}}^{l}_{l+1}$ as to the projection. However, both transfer operators ${\mathbf{I}}_{l}^{l+1}$ and ${\mathbf{P}}^{l}_{l+1}$ are technically projections.

Remark 3

As described in Section 3.2, the mesh $\pazocal{T}^{l+1}$ is obtained by uniformly refining the mesh $\pazocal{T}^{l}$. In this case, i.e. $\pazocal{T}^{l}\subset\pazocal{T}^{l+1}$, the prolongation operator ${\mathbf{I}}_{l}^{l+1}:{\pazocal{X}}^{{l}}\rightarrow{\pazocal{X}}^{{l+1}}$ obtained as described above, coincides with standard finite element interpolation.

Remark 4

In order to obtain a computationally efficient simulation, we perform the assembly of the transfer operators just once, e.g. at the beginning of the simulation. This is possible, as meshes $\{\pazocal{T}^{l}\}_{l=0}^{l=L}$ are not altered during the simulation.

6 Multilevel treatment of point-wise constraints

As discussed in Section 3.2, the minimization on level $l$ is performed with respect to some feasible set $\pazocal{L}^{l}_{c}$. The feasible set $\pazocal{L}^{l}_{c}$ is constructed as an intersection of two sets

[TABLE]

where $\pazocal{F}^{l}_{\nu}$ defines the restricted lower bound and $\pazocal{S}^{l}_{\nu}$ carries over the restricted trust region bounds. We note, that the feasible set $\pazocal{L}^{l}_{c}$ is V-cycle dependent. Therefore, it has to be re-evaluated every time the algorithm enters the given level.

In order to define the set $\pazocal{F}^{l}_{\nu}$, we follow an approach of Gelman et. al [62]. The alternatives can be found for example in [60, 54, 63]. The level dependent feasible set $\pazocal{F}^{l}_{\nu}$ is then defined as follows

[TABLE]

where the lower bound ${\mathbf{l}}^{l}_{\nu}$ is obtained component-wise as

[TABLE]

Formula (34) ensures that the prolongated coarse level corrections do not violate the fine level irreversibility constraint, thus the following relation holds

[TABLE]

Remark 5

The restriction of the point-wise constraints to the coarse level by using formula (34) is known to provide narrow approximation on the feasible set [63]. However, test runs without the coarse level constraints (34) showed similar convergence behavior, in terms of number of iterations. This can be taken as an indicator that, on the one hand, the coarse level constraints are not too restrictive and that, on the other hand, the constructed model hierarchy provides a good multi-scale decomposition of the fine level model.

The set $\pazocal{S}^{l}_{\nu}$ containing trust-region constraints inherited from the $(\nu,\mu_{1},l+1)$ iterate is also represented component-wise as

[TABLE]

where

[TABLE]

and

[TABLE]

The $\max$ and $\min$ operations in (37) - (38) are taken component-wise and ${\mathbf{e}}\in{\mathbb{R}}^{n^{l+1}}$ denotes the vector of ones. The restrictions (37) - (38) of the fine level trust region bound are weaker than the restriction of the irreversibility constraint since it is not necessary that they are satisfied by all iterates. This is due to the fact that a bounded violation of the trust region constraints is acceptable by the RMTR Algorithm 3, see [53] for the rigorous details.

Remark 6

We assume, that components of the lower bound vectors ${\mathbf{l}}^{L+1}_{\nu}$, ${\mathbf{tl}}^{L+1}_{\nu}$ equal to $-\infty$, while components of the upper bound vector ${\mathbf{tu}}^{L+1}_{\nu}$ are $\infty$.

7 Numerical experiments

In this section, we present several numerical examples, which we use in order to examine the performance of RMTR method. All of the presented examples are three dimensional and employ the length scale parameter $l_{s}=2h$, where $h$ denotes a size of the biggest element. The set of material parameters used during our simulations can be found in Table 1.

7.1 Fracture modes

At the beginning, we investigate the performance of the method on three fracture modes. We consider a three dimensional brick of size $10\times 5\times 2$ mm with $2$ mm initial notch. The specimen is then subjected to three different loading scenarios, which can force crack propagation, namely

Mode I.   - Tension (tensile stress). 2. 2.

Mode II.  - Shear (in-plane shear stress). 3. 3.

Mode III. - Tearing (out-of-plane shear stress).

Figure 5 depicts an initial set-up and details on prescribed Dirichlet boundary conditions. The final crack patterns are shown on Figure 6 by depicting the iso-surface for $c=0.95$.

7.2 Mixed mode

In our next test consider a mixed mode scenario, also know as the Nooru-Mohamed test. The numerical experiment applies both tensile and shear loads to the double-notched specimen simultaneously. Here, we utilize scheme 6a) as proposed and empirically investigated in [64]. The dimensions of the specimen, as well as loading and boundary conditions, are illustrated in the Figure 7a. In particular, we apply shear $F_{s}$ and tensile $F_{t}$ forces, such that the ratio of axial and lateral deformation remains constant throughout the test, precisely $F_{s}/F_{t}=0.5$. The Nooru-Mohamed experiment was investigated in terms of phase-field fracture models in [65].

Thanks to the specific loading protocol, the principal stresses rotate during the test, resulting in two non-interacting, non-planar, curvilinear cracks, see Figure 8. The Figure 7b illustrates the maximum principal Cauchy stress in the initial configuration.

7.3 Conchoidal fracture

The conchoidal fractures are commonly observed in fine-grained or amorphous materials such as rocks, minerals, and glasses. They are fascinating because they do not follow any natural way of the separation and result in a curved breakage, which resembles the concentric undulations of a seashell [66]. Figure 9a demonstrates an example of conchoidal fracture in quartz. From the modeling point of view, the conchoidal fracture is interesting because the crack initiates usually inside the body and not on the surface of the specimen. Therefore, the phase-field model needs to be capable of predicting the crack initiation and propagation without any initial crack or notch.

Our numerical experiment considers real rock geometry as depicted on Figure 9b. The rock mesh is unstructured and contains rough surfaces, which can only be represented well with a relatively fine mesh. Following [68], we simulate a conchoidal fracture by prescribing the displacement on the part of the upper boundary, see Figure 9b. On the same part of the boundary, the phase-field variable is enforced by prescribing Dirichlet boundary condition $c=0$. On the lower boundary, we prescribe zero Dirichlet boundary condition for both displacement and the crack phase field. The same approach, but on a parallelepipedal domain was studied by Bilgen et. al in [36]. The result of the simulation is demonstrated in Figure 10.

7.4 Pneumatic fracture

Our last numerical experiment considers pneumatic fracture. Pneumatic fracture involves the injection of air or fluid at a sufficient pressure which induces crack propagation. In order to simulate pneumatic fracturing, we follow the approach derived in [18], where a pressure term, $p$, is added into energy functional ${\Psi}({\mathbf{u}},c)$ as follows

[TABLE]

The numerical experiment was performed on a cube of size $1\times 1\times 1$ mm. Initial set-up of the simulation takes into account ten randomly distributed fractured zones, see Figure 11a. These initial cracks were specified by setting the value of the phase field to $c=1$ in those ten zones. Through the simulation, we apply zero Dirichlet boundary conditions for the displacement field on all six sides of the domain. At each pseudo-time step $t$, the pressure load is linearly increased as $p(t)=p_{0}+tp_{0}$, with an initial pressure $p_{0}=1$. The simulation time step was set-up to $\Delta t=1$s. A similar approach was considered also in [17].

Figures 11b - 11c demonstrate evolution of the fracture by illustrating crack iso-surfaces for $c=0.9$. Figure 11d depicts the obtained results by introducing a cut across the computational domain. This allows us to observe interesting crack patterns and interacting fractures inside of the computational domain.

8 Convergence and performance study

In this section, we investigate the performance of the RMTR method, Algorithm 3. Our investigation starts by analyzing the effect of different coarse level models on the performance of the RMTR. Further, we compare the RMTR method with the single level trust region solver, Algorithm 1. At the end, we examine the level-independence, complexity and scalability properties of the RMTR. While benchmarking, trust region based algorithms were set-up with the solver parameters as described in Table 2. The solution strategies terminate, when the following stoping criterion is satisfied

[TABLE]

where the current iterate is denoted by ${\mathbf{x}}_{i}$. In the context of recursive multilevel trust region algorithm, ${\mathbf{x}}_{i}$ is considered to be the fine level iterate. The criticality measure $\|\pazocal{E}({\mathbf{x}}_{i},f)\|$ is as defined in (49) and $\|{\mathbf{x}}_{i-1}-{\mathbf{x}}_{i}\|$ denotes the correction size. The RMTR solver was configured as a V-cycle with one pre/post-smoothing step and two coarse solves. We solved the trust region subproblems on levels $l=1,\dots,L$ approximately with 10 steps of the projected Conjugate-Gradient method preconditioned with Jacobi preconditioner [25]. On the coarsest level, we employed the active set strategy [69] with the sparse linear solver MUMPS [70].

Our implementation of the parallel multilevel phase-field fracture simulations is based on the finite element framework MOOSE [46]. Presented solution strategies were implemented as part of our open-source, C++ library Utopia [71]. Utopia is an embedded domain specific language (EDSL), designed to make implementations for parallel computing as transparent as possible. It enables us to write complex scientific code by using expressions similar to Matlab, whereas the complexity of parallelization is hidden in different back-ends. The presented results were computed with a PETSc [72] back-end.

All experiments were performed on the local cluster at the Institute of Computational Science (ICS), Università della Svizzera Italiana. The cluster consists of $42$ compute nodes, each equipped with 2 Intel R E5-2650 v3 processor with a clock frequency of 2.60GHz. Having 10 cores per processor leads to 20 cores per node. Memory per node is 64 GB.

8.1 Effect of different coarse level models on the performance of the RMTR

In this section, we investigate the performance of the RMTR method configured with different coarse level models, as described in Section 4. We monitor number of nonlinear V-cycles required by RMTR method to converge to the desired tolerance during the time step, see Figure 12. As we can see from Figure 12, all variants of RMTR method converge to the desired tolerance. However, convergence speed differs. Based on our numerical experiments, we conclude that the 1st order consistency approach (18) yields the worst convergence rates. This is especially true in the post-stage, where fracture occurs and propagates. As mentioned in Section (4.1), the slow convergence is connected to the fact that the coarse level models are not suitable representations of the fine level fracture zones.

Figure 12 demonstrates that all three approaches: Galerkin, 1st order consistency and 2nd order consistency achieve worse convergence rates as our novel solution dependent 2nd order model (23). We can compute achieved speed-up with respect to any other coarse level model as follows

[TABLE]

From the obtained results, Table 3, we can clearly see that usage of our solution dependent 2nd order coarse level model is beneficial and leads to significant speed-up. For example, for the mixed mode test, the speed-up equals $7.42$, when compared to 1st order consistency approach.

All results reported in the reminder of this section employ our RMTR method with the coarse level models based on the solution dependent second order consistency approach, Section 4.4.

8.2 Single level TR vs. RMTR

In this section, we compare the performance of RMTR Algorithm 3 with its single level variant, TR Algorithm 1. To this aim, we measure the number of nonlinear iterations required by TR method vs. the number of nonlinear V-cycles required by RMTR method. Figures 13 - 14 depict obtained results for different time-steps. As we can observe from Figures 13 - 14, the number of nonlinear iterations grows once the simulation reaches the post-stage phase, where the crack finally propagates. This behavior is not surprising, as in the post-stage phase the non-convexity and the non-linearity become stronger. For example, from the tension test in Figure 13a, we can see that crack propagation starts around time $t=4.5\times 10^{-3}$.

Our numerical experiences indicate that the RMTR method always outperforms the single level TR method. On the average, we observe a reduction in the number of required iterations approximately by factor of ten. The biggest decrease can be observed for the pressure example, where the number of iterations required by RMTR method is 32 times lower than the number of iterations required by the single level TR method. Another interesting results were obtained for the mixed mode example, see Figure 14a. The single level TR algorithm encounters a failure, as the solver did not converge within $10,000$ iterations. At the same time, the multilevel variant converges without any problems. However, we have to remark that the single level TR method converges, while performing the simulation with smaller time-step.

8.3 Level-independence

The classical multilevel solvers, for instance multigrid are level-independent, in the sense, that the number of iterations stays constant with an increased number of dofs. Here, we investigate the level-independence of our RMTR algorithm by measuring number of required V-cycles for three fracture modes while increasing number of levels/dofs. The Table 4 provides number of dofs used on each level of multilevel hierarchy as well as respective length-scale parameter $l_{s}$. It is of major importance to notice that the length-scale parameter $l_{s}$ decreases with increased resolution. As a consequence, volumetric approximation of the fracture reassembles sharp fracture surface more closely as refinement level increases.

Figure 15 shows, that the number of required V-cycles increases slightly with the number of levels. This might be caused by two factors: Firstly, phase-field fracture problems discretized on different grid are not able to represent the same fractures. Therefore, simulating fracture propagation on different discretization level does not imply that we are solving the same problem. The second reason, why we are not able to observe level-independence might come directly from the use of the trust region globalization. During certain stages of the simulation, the trust region radius might become too small. In this case, the sum of corrections coming from upper levels of the multilevel hierarchy already satisfy the constraints provided by the trust region radius on the finest level. As a consequence, recursion terminates before reaching the coarsest level. The low frequencies of the error are not eliminated, which results in the higher number of nonlinear V-cycles. However, this is a very small price to pay in order to obtain a globally convergent multilevel solution strategy.

In comparison, Figure 16 demonstrates how number of required nonlinear iterations changes with increased refinement level for the single level TR method. We note, that growth in number of iterations is prevalent from one refinement level to the another. For example, for the tension example, Figure 16a, the accumulated number of iterations increases by factor $3.17$ from refinement level $5$ to refinement level $6$.

8.4 Complexity and scalability

One of the most appealing property of any the multilevel algorithm is its optimal complexity $\mathscr{O}(n)$. This means, that the required solution time increases linearly with increased number of dofs. We can demonstrate the computational complexity by measuring the execution time of one V-cycle while increasing number of dofs. Our experiment considers tension test from 7.1. The obtained results, depicted on Figure 17a, demonstrate that our implementation of the RMTR method is indeed of the optimal complexity. One should note, that the reported times are based on the results obtained on a serial machine.

Furthermore, we examine the strong scalability properties of our Utopia implementation of the RMTR method. We consider tension test with fixed number of dofs, namely $782,340$. Then, we measure required solution time while increasing number of cores. This allows us to compute the relative speedup by using following formula

[TABLE]

where $T_{20}$ represents computational time required by 20 cores and $T_{p}$ symbolizes the computational time required by $p$ cores. Obtained results are reported on Figure 17b, which demonstrates that our implementation gives rise to almost ideal scaling up to 250 cores. For more than $250$ processors, scaling is not ideal anymore, which is caused by an insufficient amount of dofs on the coarsest level. Further examination of scaling properties of the RMTR method as well as performing tests with larger examples (with more dofs) is left for the future work.

Remark 7

The measured time for both complexity and scalability tests contains also time spent for finite-element assembly.

9 Conclusion

In this work, we proposed a non-linear multi-level method for solving the monolithic systems arising from the phase-field fracture simulations. Our method is based on recursive multi-level trust-region methods (RMTR), see [37, 38]. The proposed variant of RMTR method employs novel level dependent objective functions, which combine coarse level discretizations with a fine level representation of the fracture. The usage of this novel formulation improved the performance of the RMTR method substantially. Furthermore, a comparison with single level trust region method was made and demonstrated a reduction in the number of iterations approximately by order of magnitude. Performed experiments also indicated that the our implementation of the RMTR method is of optimal complexity and shows good scaling behavior.

Appendix A Euler-Lagrange equations

The first order necessary conditions for the minimization problem (5) are defined as:

Find pair $({\mathbf{u}},c)\in\bm{V}^{t}\times H^{1}(\Omega)$, such that

[TABLE]

where $\bm{V}^{t}:=\{{\mathbf{u}}\in\bm{H}^{1}(\Omega)\ |\ {\mathbf{u}}=\bm{g}^{t}\ \text{on}\ \ \partial\Omega_{D}\}$ and $W^{t}:=\{c\in H^{1}(\Omega)\ |\ 0\leq c^{t-1}\leq c\leq 1\ \text{a.e.\ on}\ \Omega\}$. The symbol $c^{t-1}$ denotes the solution of the phase-field on the previous time-step. The directional derivatives of $\Psi({\mathbf{u}},c)$ with respect to ${\mathbf{u}}$ and $c$, denoted by $\nabla_{{\mathbf{u}}}\Psi$ and $\nabla_{c}\Psi$, respectively are defined as

[TABLE]

where $k\approx 0$ denotes a regularization parameter.

Appendix B Additive decomposition of elastic energy

Following [6], we additively decompose the elastic energy density $\psi_{e}(\bm{\varepsilon})$ into $\psi_{e}^{+}(\bm{\varepsilon})$ due to the tension and $\psi_{e}^{-}(\bm{\varepsilon})$ due to the compression by using following relation

[TABLE]

The symbols $\mu$ and $\lambda$ in (45) denote the Lamé constants. The bracket symbol is defined as \big{<}x\big{>}^{\pm}:=\frac{(|x|\pm x)}{2}. The strain tensor related to the desirable mode (tension, or compression) is obtained by the spectral decomposition of the principal strains

[TABLE]

where $\bm{n}_{i}$ are the principal directions. The Cauchy stress tensor $\boldsymbol{\sigma}:=\frac{\partial\psi}{\partial{\bm{\varepsilon}}}$ is also decomposed in an additive manner as

[TABLE]

with

[TABLE]

Appendix C Projection onto feasible set

Let us consider feasible set $\pazocal{F}:=\{{\mathbf{x}}\in{\mathbb{R}}^{n}\ |\ {\mathbf{l}}\leq{\mathbf{x}}\leq{\mathbf{u}}\}$, where the inequalities are meant component-wise. The projection $\mathscr{P}({\mathbf{y}})$ of vector ${\mathbf{y}}$ onto $\pazocal{F}$ is also defined component-wise as

[TABLE]

for $\forall k\in\{1,\dots,n\}$.

Acknowledgments

The authors thank Maria Nestola, Patrick Zulian and Hardik Kothari for many useful discussions and suggestions. The authors thank Cyrill von Planta for providing rock geometry. This work was supported by Swiss National Science Foundation (SNF) and DFG (Deutsche Forschungsgemeinschaft) under the project SPP1748 "Reliable Simulation Techniques in Solid Mechanics. Development of Non-standard Discretization Methods, Mechanical and Mathematical Analysis".