Key words. adaptive finite element method, brittle fracture, free-discontinuity ... to J(f), f± are the inner and outer traces of f on J(u) with respect to νf , and Dcf is ...
AN ADAPTIVE FINITE ELEMENT APPROXIMATION OF A VARIATIONAL MODEL OF BRITTLE FRACTURE ¨ SIOBHAN BURKE , CHRISTOPH ORTNER , AND ENDRE SULI Abstract. The energy of the Francfort–Marigo model of brittle fracture can be approximated, in the sense of Γ-convergence, by the Ambrosio–Tortorelli functional. In this work we formulate and analyze two adaptive finite element algorithms for the computation of its (local) minimizers. For each algorithm we combine a Newton-type method with residual-driven adaptive mesh refinement. We present two theoretical results which demonstrate convergence of our algorithms to local minimizers of the Ambrosio–Tortorelli functional. Key words. adaptive finite element method, brittle fracture, free-discontinuity problem, Ambrosio–Tortorelli functional AMS subject classifications. 65N30, 65N50, 74R10
1. Introduction. Beginning with the work of Francfort and Marigo [24], the mathematical theory of quasi-static brittle fracture mechanics has experienced a rapid and successful development. Upon recasting Griffith’s idea of balancing energy release rate with a fictitious surface energy [26] as an energy minimization problem, Francfort and Marigo were able to formulate a model that was free of the usual constraints of fracture mechanics such as a predefined and piecewise smooth crack path. With the help of the theory of free-discontinuity problems [5], this model was soon shown to be well-posed in a surprisingly general setting [17, 18, 23]. We briefly review the model in Section 1.1. The model of Francfort and Marigo is posed in terms of the minimization of a highly irregular energy functional, which also occurs in image segmentation where it is known as the Mumford– Shah functional [29]. Several methods have been proposed in the literature, which regularize this energy in order to render the problem accessible to numerical simulation [10, 15, 30]. Such methods typically use the theory of Γ-convergence [12] to construct approximating functionals whose minimizers converge to those of the original functional. In our experience, the Ambrosio–Tortorelli approximation [3, 4] is one of the most promising approaches. A particularly nice feature of the Ambrosio–Tortorelli functional is that its minimization can be reduced to the solution of elliptic boundary value problems that are straightforward to discretize, for example, by a finite element method. This approach has been used successfully by Bourdin et al. [8, 9, 11] for the simulation of problems that are usually inaccessible to traditional methods. A brief review of the Ambrosio–Tortorelli approximation is given in Section 1.2. The Ambrosio–Tortorelli approximation can be understood as a phase field model for the crack set. To resolve the phase field variable, the mesh near the crack has to be significantly finer than the mesh that would be required to resolve the elastic deformation. Since we do not know the crack path in advance, it is therefore a natural idea to use an adaptively refined mesh. Using established techniques of adaptive finite element theory [35], we derive a residual estimate for the gradient of the Ambrosio–Tortorelli functional (see Section 3.1). In Section 3.2, we formulate two adaptive algorithms for which we use the residual estimates to steer the mesh refinement. In Section 4 we prove that the second of these algorithms converges to a critical point of the Ambrosio–Tortorelli functional. We also prove that the first algorithm converges to a critical point under the assumption of the associated residual converging to zero. We conclude with an implementation of the adaptive algorithms for a computational example in Section 5. In order to lay out the main ideas, our analysis in the present work is restricted to linearized elasticity in anti-plane deformation, and to linear finite elements. We will extend our results to more general approximations and a wider range of models in future work. 1.1. The Francfort–Marigo Model of Brittle Fracture. In order to introduce the Francfort–Marigo model of brittle fracture, we briefly define the space of special functions of bounded variation [5, 21]. A knowledge of this space is not necessary for the main ideas contained in the paper. 1
For p ∈ [1, ∞] and Ω an open domain in RN , we use Lp (Ω) to denote the standard Lp -spaces and H1 (Ω) to denote the standard Hilbertian Sobolev space. The N -dimensional Lebesgue and Hausdorff measures are denoted by LN and HN respectively. We say that a function f ∈ L1 (Ω) has bounded variation, if nZ o sup f divϕ dx : ϕ ∈ C10 (Ω; RN ), |ϕ| ≤ 1 < ∞. Ω
The space of functions of bounded variation is denoted BV(Ω). A function of bounded variation can exhibit discontinuities that are reflected in its distributional gradient. Given a function f ∈ BV(Ω) the distributional derivative of f , denoted Df , is a Radon measure, which can be decomposed as Df = ∇f LN + (f + (x) − f − (x)) ⊗ νf (x)HN −1 ⌊J(f ) + Dc f, where ∇f is called the approximate gradient of f , J(f ) is the jump set of f , νf is the unit normal to J(f ), f ± are the inner and outer traces of f on J(u) with respect to νf , and Dc f is called the Cantor part of the derivative. We refer to [5] for the precise definitions of these terms. The space of special functions of bounded variation is then defined as SBV(Ω) := {f ∈ BV(Ω) : Dc f = 0}. We are now in a position to describe the Francfort–Marigo model of brittle fracture [24]. The crack-free reference configuration of a linearly elastic body is denoted by Ω. The set Ω is taken to be an open, bounded and connected domain with Lipschitz boundary ∂Ω (we shall later on relax this assumption). For each u ∈ SBV(Ω), and for each Hausdorff measurable set Γ, we define the bulk, surface, and total energy, respectively, by Z |∇u|2 dx, (1.1) EB (u) := Ω\J(u)
ES (Γ) := κHN −1 (Γ), and EB (u) + ES (Γ), E(u, Γ) := +∞,
(1.2) if HN −1 (J(u) \ Γ) = 0, otherwise.
(1.3)
The energy functional E(u, Γ) reflects Griffith’s principle that, to create a crack one has to spend an amount of elastic energy that is proportional to the area of the crack created [26]. The constant of proportionality κ is called fracture toughness and is often also denoted Gc (critical energy release rate). The crack set Γ and the jump set J(u) are decoupled in the definition of the total energy in order to be able to impose irreversibility of the crack evolution. The quasistatic evolution will be obtained upon minimizing E subject to several constraints, which we introduce momentarily. We wish to study how the body evolves in time under the action of a varying load g(t), which is applied on an open subset ΩD ⊂ Ω of positive N -dimensional Lebesgue measure. We assume that g ∈ L∞ (0, T ; W1,∞ (Ω)) ∩ W1,1 (0, T ; H1 (Ω)), and we define A(g(t)) := u ∈ SBV(Ω) : u|ΩD = g(t)|ΩD .
The fact that the Dirichlet condition is imposed on a set of positive N -dimensional Lebesgue measure is mostly technical and ensures that the jump set on the Dirichlet boundary ∂ΩD ∩ Ω is well-defined. Since we work in anti-plane deformation the boundary displacement is applied in a direction perpendicular to the plane in which the initial configuration of the domain lies. We call ∂ΩN := ∂Ω\∂ΩD the Neumann boundary. It is most intuitive to introduce the Francfort–Marigo model through a time discretization. Let 0 = t0 < t1 < · · · < tF = T be a discretization of the time interval [0, T ], with ∆t := max {tk − tk−1 : k = 1, . . . , F }. Given an initial crack Γ(0) (which should be the jump set of an SBV function), we seek (u(tk ), Γ(tk )), k = 1, . . . , F , such that 2
1. Irreversibility: Γ(tk ) ⊃ Γ(tk−1 ); 2. Global stability: E(u(tk ), Γ(tk )) ≤ E(v, Γ) ∀v ∈ A(g(tk )) and Γ ⊃ Γ(tk−1 ). In practise, this formulation requires the successive solution of the global minimization problems u(tk ) := argmin EB (v) + ES (J(v)) ∪ Γ(tk−1 )), v∈A(g(tk ))
Γ(tk ) := J(u(tk )) ∪ Γ(tk−1 ).
(1.4)
In this formulation the problem is accessible to the direct method of the calculus of variations to prove the existence of solutions to the time-discrete Francfort–Marigo model [1, 2, 5, 19]. It remains to describe the limit of the time-discrete evolution as ∆t → 0. In full generality, it was first shown by Francfort and Larsen [23] (an earlier result is due to Dal Maso and Toader [18] and the extension to finite elasticity was developed by Dal Maso, Francfort, and Toader [17]) that it is indeed possible to extract weakly convergent subsequences of the discrete solutions and thus prove the existence of a trajectory u ∈ BV(0, T ; SBV(Ω)) with crack set Γ(t) = ∪s0 satisfies the following Γ-convergence result [11]. Let us define Iε (ˆ u, vˆ), if uˆ ∈ H1 (Ω), uˆ = g(t) on ΩD , vˆ ∈ H1 (Ω; [0, 1]), Gε (ˆ u, vˆ) := +∞, otherwise, E(ˆ u, J(ˆ u)), if uˆ ∈ A(g(t)), vˆ = 1 a.e. on Ω, G(ˆ u, vˆ) := +∞, otherwise; then Gε Γ-converges to G in L1 (Ω) × L1 (Ω) as ε → 0. 3
This is a modification of the original Ambrosio–Tortorelli Γ-convergence result [3] for the approximation of the Mumford–Shah functional [29]. The existence of minimizers to Iε has been shown in [4] for each ε, η > 0. Further generalisations of the Γ-convergence result have since been established in [7, 22]. Using the Ambrosio–Tortorelli functional, an approximation to the time-discrete Francfort– Marigo model can be computed as follows. At time t = t0 , find (uε (t0 ), vε (t0 )) ∈ argmin{Iε (ˆ u, vˆ) : uˆ ∈ H1 (Ω), u ˆ = g(t0 ) on ΩD ; vˆ ∈ H1 (Ω; [0, 1])}. At subsequent times t = tk , k = 1, . . . , F , find (uε (tk ), vε (tk )) satisfying (uε (tk ), vε (tk )) ∈ argmin{Iε (ˆ u, vˆ) : u ˆ ∈ H1 (Ω), u ˆ = g(tk ) on ΩD ; vˆ ∈ H1 (Ω), vˆ ≤ vε (tk−1 )}.(1.6) At the fixed time t = tk , the function uε (tk ) is an approximation of the displacement of the body u(tk ), with uε (tk ) → u(tk ) in L1 (Ω) as ε → 0. The auxiliary function vε (tk ) is introduced in order to represent the crack Γ(tk ). A truncation argument shows that 0 ≤ vε (x, tk ) ≤ 1 a.e. x ∈ Ω. The crack is then approximated by the subset of the domain on which vε (tk ) takes values close to zero. Conversely, the unfractured part of the body is represented by the subset of the domain on which vε (tk ) takes values close to one. The transition layer between the two regimes has thickness of order ε. In the limit ε → 0 the minimization of Iε requires that vε (tk ) → 1 almost everywhere. Consequently, the transition layer becomes infinitely thin, with vε (tk ) → 0 only on the (N − 1)-dimensional surface representing the crack. It was shown by Giacomini [25] that an evolution satisfying (1.6) converges as ∆t, ε → 0 to an evolution satisfying the conditions of the time-continuous Francfort–Marigo model. We will therefore restrict our consideration to the problem of minimizing the Ambrosio–Tortorelli functional at the fixed moment in time t = tk , and for fixed values of ε and η. The condition vˆ ≤ vε (tk−1 ) enforces the irreversibility of the crack [25]. In practise, however, we choose to implement the irreversibility criterion through the following equality constraint introduced by Bourdin [8]. If at time t = tk , k ∈ {0, 1, . . . , F − 1}, the set ¯ : vε (x, tk ) < CRTOL} CR(tk ) := {x ∈ Ω
(1.7)
is non-empty for some small specified tolerance CRTOL, then fix vε (x, ti ) = 0 for all x ∈ CR(tk ) and all i such that k < i ≤ F. Thus, if at a particular time t = tk , vε (x, tk ) is close enough to zero to indicate that the point x lies on the crack path, then vε is set to be zero at that point for all subsequent time steps. This simplifies the minimization over vˆ by allowing the use of an unconstrained minimization algorithm. However, this modification of the irreversibility condition vˆ ≤ vε (tk−1 ) is not as yet shown to be equivalent to irreversibility of the crack as ∆t, ε → 0. A numerical discretization of this minimization problem has been proposed and implemented by Bourdin et al. [8, 9, 11] in which a finite element method is used in conjunction with a minimization algorithm to compute minimizers of the Ambrosio–Tortorelli functional. The method uses a fine mesh to discretize the whole domain. It does, however, seem apparent to us that the behaviour of minimizers can be sufficiently resolved using a mesh which is fine only in a layer around the crack, whilst remaining relatively coarse elsewhere. Naturally, we do not know the location of the crack a priori. Thus, we propose to use an adaptive finite element method for the numerical computation of minimizers. One could argue that the refinement of only part of the mesh causes the crack to favour growth in this direction; we believe, however, that a carefully designed adaptive algorithm can circumvent this eventuality. We refer to Remark 4 for a further discussion of this point. 4
1.3. Critical Points. In the Ambrosio–Tortorelli model, Lipschitz regularity of the domain is not required. Since, in practise, it is more convenient to model a pre-existing crack by a slit domain rather than an initial crack field v0 we will, from now on, drop the assumption that Ω is a Lipschitz domain. Motivated by the fact, that we will need to partition Ω for the purpose of defining a finite element approximation, we shall assume that Ω is a polyhedral domain. By this, we simply mean that Ω possesses a finite partition into non-degenerate N -simplices: there exist open, non-degenerate simplices T1 , . . . , TK ⊂ Ω such that LN (Ω \ ∪k Tk ) = 0 (see also Section 2). In that case, it is clear that the usual trace and embedding theorems for Sobolev spaces hold. Since we consider the minimization of the Ambrosio–Tortorelli functional at the fixed time t = tk , it is useful to define the following function spaces: H1D (Ω) := {w ∈ H1 (Ω) : w = g(tk ) on ΩD }, H1D0 (Ω) := {w ∈ H1 (Ω) : w = 0 on ΩD }.
(1.8) (1.9)
We also fix ε and η, and for simplicity set κ = 1. Accordingly we choose to relabel the Ambrosio– Tortorelli functional as I : H1D (Ω) × H1 (Ω) → R ∪ {+∞} where Z 2 (1.10) (v + η)|∇u|2 + α(1 − v)2 + ε|∇v|2 dx, I(u, v) := Ω
and α = 1/4ε. It can be seen, using a truncation argument, that any local minimizer (u, v) of I (in the H1 (Ω) × H1 (Ω) topology) satisfies 0 ≤ v(x) ≤ 1 a.e. in Ω. Thus all relevant test functions for v lie in the space L∞ (Ω). As such, in the following discussion of differentiability of I, we work with test functions for v from the space H1 (Ω) ∩ L∞ (Ω). Proposition 1.1. I is Fr´echet-differentiable in H1 (Ω) × (H1 (Ω) ∩ L∞ (Ω)). Proof. Let (u, v) ∈ H1 (Ω) × (H1 (Ω) ∩ L∞ (Ω)). A straight forward calculation yields that the directional derivative of I at (u, v) in the direction (ϕ, ψ) ∈ H1 (Ω) × (H1 (Ω) ∩ L∞ (Ω)) is given by Z Z vψ|∇u|2 + α(v − 1)ψ + ε∇v · ∇ψ dx, I ′ (u, v; ϕ, ψ) = 2 (v 2 + η)∇u · ∇ϕ dx + 2 Ω
Ω
=: 2a(v; u, ϕ) + 2b(u; v, ψ).
This is in fact the Fr´echet derivative since g(u, v; ϕ, ψ) := I(u + ϕ, v + ψ) − I(u, v) − I ′ (u, v; ϕ, ψ) Z = [ψ 2 |∇u|2 + |∇ϕ|2 (v + ψ)2 + η + (4vψ + 2ψ 2 )∇u · ∇ϕ + αψ 2 + ε|∇ψ|2 ] dx, Ω
is such that |g(u, v; ϕ, ψ)| ≤ kψk2L∞ (Ω) k∇uk2L2 (Ω) + k(v + ψ)2 + ηkL∞ (Ω) k∇ϕk2L2 (Ω)
+ 4kvkL∞(Ω) kψkL∞ (Ω) k∇ukL2 (Ω) k∇ϕkL2 (Ω) + 2kψk2L∞ (Ω) k∇ukL2 (Ω) k∇ϕkL2 (Ω)
+ αkψk2L2 (Ω) + εk∇ψk2L2 (Ω) , and thus
|g(u, v; ϕ, ψ)| → 0 as (kϕkH1 (Ω) + kψkH1 (Ω) + kψkL∞ (Ω) ) → 0. kϕkH1 (Ω) + kψkH1 (Ω) + kψkL∞ (Ω) Remark 1. We note that I(u, v) is not finite for all (u, v) ∈ H1 (Ω) × H1 (Ω), and thus I is not Gateaux-differentiable in H1 (Ω) × H1 (Ω). This motivates the following definition of a critical point. Definition 1.2. We say that (u, v) ∈ H1D (Ω) × (H1 (Ω) ∩ L∞ (Ω)) is a critical point of I if ′ I (u, v; ϕ, ψ) = 0 for all ϕ ∈ H1D0 (Ω) and ψ ∈ H1 (Ω) ∩ L∞ (Ω). 5
Proposition 1.3. If (u, v) ∈ H1D (Ω) × (H1 (Ω) ∩ L∞ (Ω)) is a critical point of I then 0 ≤ v(x) ≤ 1 for a.e. x ∈ Ω. Proof. Suppose that (u, v) is a critical point of I and that there exist subsets A1 and A2 of Ω with |A1 ∪ A2 | > 0 such that v > 1 on A1 and v < 0 on A2 . Since (u, v) is a critical point it follows that Z vψ|∇u|2 + α(v − 1)ψ + ε∇v · ∇ψ dx = 0 ∀ψ ∈ (H1 (Ω) ∩ L∞ (Ω)). Ω
Choosing
yields Z
A1
1 − v(x) 0, ψ(x) = −v(x),
2
2
2
v(1 − v)|∇u| − α(1 − v) − ε|∇v|
if x ∈ A1 , if x ∈ Ω\(A1 ∪ A2 ), if x ∈ A2 ,
dx −
Z
A2
2 v |∇u|2 + αv(v − 1) + ε|∇v|2 dx = 0,
which is a contradiction since the left-hand side is strictly negative. 2. Finite Element Discretization. 2.1. Finite Element Discretization. Since we assumed that Ω is a polyhedral domain (see Section 1.3), we may discretize it as follows. Let Th be a subdivision of Ω into N -dimensional ¯ = ∪T ∈T T¯ and Ti ∩ Tj = ∅ for Ti , Tj ∈ Th with i 6= j. The open simplexes T ∈ Th , such that Ω h subdivision Th is chosen in such a way that the boundary of ΩD is discretized as the union of faces of simplices from Th . We define h := maxT ∈Th diam(T ) and each simplex T ∈ Th is taken to be an affine transformation of the open unit simplex T¯ := {¯ x = (¯ x1 , . . . , x ¯N ) : 0 < x ¯i , i = 1, . . . , N, 0 < x ¯1 + · · · + x ¯N < 1}. Each simplex T ∈ Th is called an element. We assume that the subdivision is conforming, that is, the intersection of the closure of any two elements is either empty or is along an entire kdimensional face, 0 ≤ k ≤ N − 1. We also require that the subdivision is shape-regular, i.e., sup T ∈Th
hT ≤ K, ρT
for some K ∈ (0, ∞), where hT := diam(T ) and ρT is the diameter of the largest N -dimensional ball contained in T . Let Nh ⊂ N denote the set of vertices in the subdivision. For each vertex i ∈ Nh , let xi denote the position of the ith vertex and let ζi be the continuous piecewise linear basis function such that ¯ D }. ζi (xj ) = δij . Define ND,h := {i ∈ Nh : xi ∈ Ω Let Eh denote the set of (N − 1)-dimensional faces in the subdivision with ED,h := {e ∈ Eh : ¯ D }, EN,h := {e ∈ Eh : e ⊆ ∂ΩN } and EI,h := {e ∈ Eh \(ED,h ∪ EN,h )}. We define ED,h , EN,h e⊆Ω and EI,h as the union of all faces in ED,h , EN,h and EI,h respectively. Let ωx be the union of elements T ∈ Th that have x as the position of a vertex. For a face e ∈ Eh and element T ∈ Th define ωe := ∪x∈e ωx , ωT := ∪x∈T ωx and he = diam(e). We now define the following finite element spaces: n X o Xh = λi ζi : λi ∈ R , i∈Nh
Xh,D0 =
n X
i∈Nh
o λi ζi : λi ∈ R, λi = 0 ∀i ∈ ND,h , 6
and the finite element space at time t = tk : n X o k λi ζi : λi ∈ R, λi = g(tk , xi ) ∀i ∈ ND,h . Xh,D = i∈Nh
k Since we consider a fixed time tk we relabel Xh,D as Xh,D for ease of notation. Also, for simplicity, we assume throughout that g lies in the finite element space Xh . For our subsequent analysis, we require the discrete formulation to satisfy a maximum principle analogous to Proposition 1.3. In order to accomplish this we use a mass lumping approximation [34, Chapter 11] for I together with an assumption on the stiffness matrix, which can be achieved through typical mesh regularity conditions. The mass lumping approximation of I is defined to be the functional: Z (2.1) Ph (vh2 ) + η |∇uh |2 + αPh (vh − 1)2 + ε|∇vh |2 dx, Ih (uh , vh ) := Ω
where Ph is the standard nodal interpolation operator [13, Section 3.3]. Let K be the |Nh | × |Nh | stiffness matrix with entries (kij )i,j∈Nh , then we assume that Z ∇ζi · ∇ζj dx ≤ 0 ∀i 6= j ∈ Nh . (2.2) kij := Ω
This condition has been studied in detail in two-dimensions [16], [33, p.78] and three-dimensions [27]. We note that this assumption is only made for technical purposes and in practise we do not believe that it is necessary for Proposition 2.2 to hold. Definition 2.1. We say that (uh , vh ) ∈ Xh,D ×Xh is a critical point of Ih if, Ih′ (uh , vh ; ϕh , ψh ) = 0 for all ϕh ∈ Xh,D0 and ψh ∈ Xh , where Ih′ (uh , vh ; ϕh , ψh ) = 2ah (vh ; uh , ϕh ) + 2bh (uh ; vh , ψh ), Z ah (vh ; uh , ϕh ) := Ph (vh2 ) + η ∇uh · ∇ϕh dx, ZΩ Ph (vh ψh )|∇uh |2 + αPh (vh − 1)ψh ) + ε∇vh · ∇ψh dx. bh (uh ; vh , ψh ) := Ω
Proposition 2.2. Suppose that (uh , vh ) ∈ XhD × Xh such that bh (uh ; vh , ψh ) = 0 for all ψh ∈ Xh then 0 ≤ vh (x) ≤ 1 for all x ∈ Ω. Proof. Let (uh , vh ) ∈ XhD × Xh be such that bh (uh ; vh , ψh ) = 0 for all ψh ∈ Xh . Let P vh = j∈Nh vj ζj and suppose, for contradiction, that there exist subsets J1 and J2 of Nh with vj > 1 for all j ∈ J1 and vj < 0 for all j ∈ J2 . First suppose that J1 is non-empty and i ∈ J1 is such that vi ≥ vj for all j ∈ Nh . Consider the patch of elements Bi := supp(ζi ), and define Mi := {j ∈ Nh : j ∈ Bi }. On taking ψh = ζi we have bh (uh ; vh , ζi ) = 0, which implies that Z Z Z 2 Ph ((vh − 1)ζi ) dx Ph (vh ζi )|∇uh | dx − α ∇vh · ∇ζi dx = − ε Bi
Bi
Bi
vi X α =− |∇uh |2T |T | − (vi − 1)|Bi | 3 3 T ⊂Bi
< 0. Since vh = Z
P
Bi
j∈Nh
(2.3)
vj ζj , on noting that the rows of the stiffness matrix sum to zero, we have
∇vh · ∇ζi dx =
X
j∈Mi
kij vj =
X
j∈Mi
kij (vj − vi ) +
X
kij vi =
j∈Mi
Thus, using the assumption that kij ≤ 0 for i 6= j, it follows that contradicts (2.3). 7
X
j∈Mi
R
Bi
kij (vj − vi ).
∇vh · ∇ζi dx ≥ 0, which
A similar argument can be used to contradict the existence of vertices in J2 . In an analogous way to the definition of CR(tk ), for a given vh ∈ Xh and a small specified tolerance CRTOL we define CRh (tk ) := {j ∈ Nh : vh (xj , tk−1 ) < CRTOL}. In the finite element approximation, at time t = tk , we seek to find (uh , vh ) ∈ argmin{Ih (ˆ uh , vˆh ) : u ˆh ∈ Xh,D , vˆh ∈ Xh with vˆh (xj ) = 0 for all j ∈ CRh (tk )}.
(2.4)
2.2. The Alternate Minimization Algorithm. The minimization of the functional Ih is a nontrivial task since the term Ph (vh2 ) |∇uh |2 renders the functional nonconvex. A number of minimization schemes can be employed; we note however that none would in general be able to find the location of global minimizers, at least not easily. Instead we must be satisfied with being able to locate local minimizers. (In any case, as we have previously noted, it is unclear to us that finding global minimizers of the Francfort–Marigo model is physically justified). The minimization will be achieved using an alternate minimization algorithm proposed by Bourdin et al. [11]. We state the algorithm for the minimization of I over the infinite-dimensional space H1D (Ω)×H1 (Ω) at time t = tk . The idea is as follows. Although the functional I is nonconvex with respect to the pair (u, v), I is convex with respect to u and v separately. Thus it is a straightforward computation to minimize with respect to one variable at a time. For some termination tolerance VTOL, the algorithm proceeds as follows: 1. Set v 0 = 1 if k = 0 and v 0 = v(tk−1 ) if k > 0. 2. For n = 1, 2, . . . • un = argmin {I(z, v n−1 ) : z ∈ H1D (Ω)} • v n = argmin {I(un , w) : w ∈ H1 (Ω), w = 0 on CR(tk−1 )} • Repeat until kv n − v n−1 kL∞ (Ω) < VTOL 3. Set u(tk ) = un and v(tk ) = v n . The algorithm has been successfully implemented by Bourdin et al. [8,9,11] for a range of computational examples. However, to date there does not exist a proof of convergence for the algorithm to a minimizer as n → ∞. It appears to us that the proof of [8, Theorem 1] is incomplete, even thought the result is certainly correct. The missing arguments can be found in the present paper, in the proofs of Theorems 4.1 and 4.2 (in particular step 4 in the proof of Theorem 4.1 and steps 3 and 4 in the proof of Theorem 4.2). For a proof of local convergence to isolated local minimizers see [8, Theorem 2]. Remark 2. The alternating minimization algorithm can be understood as a Newton-type descent method with a special starting guess. To see this, let u0 ∈ H1D (Ω) and let v 0 = argminH1 (Ω) I(u0 , ·). The first step of the alternating minimization algorithm, computing the pair (u1 , v 1 ), can be written as ∂u I(u1 , v 0 ) = 0
and
∂v I(u1 , v 1 ) = 0,
(2.5)
where ∂u I and ∂v I denote the partial derivatives of I. Since I is quadratic in each coordinate, with non-singular partial derivatives ∂uu I and ∂vv I, the equalities (2.5) are equivalent to ∂uu I(u0 , v 0 )(u1 − u0 ) = −∂u I(u0 , v 0 )
and
∂vv I(u1 , v 0 )(v 1 − v 0 ) = −∂v I(u1 , v 0 ).
It can now be easily seen that these two equations are equivalently written as 1/2 ∂uu I(u0 , v 0 ) 0 u ¯ − u0 ∂u I(u0 , v 0 ) = − and 0 ∂vv I(u0 , v 0 ) ∂v I(u0 , v 0 ) v¯1/2 − v 0 1 ∂uu I(¯ u1/2 , v¯1/2 ) 0 u ¯ −u ¯1/2 ∂u I(¯ u1/2 , v¯1/2 ) = − , 0 ∂vv I(¯ u1/2 , v¯1/2 ) v¯1 − v¯1/2 ∂v I(¯ u1/2 , v¯1/2 ) where (¯ u1/2 , v¯1/2 ) := (u1 , v 0 ) and (¯ u1 , v¯1 ) := (u1 , v 1 ). 8
3. Adaptive Algorithm. The nature of solutions to the minimization problem strongly motivates the use of an adaptive finite element method for their computation. Such methods use a local refinement indicator, based on the computed solution, to identify those elements where mesh refinement would be most beneficial for improving the accuracy of the solution. It is now a well-established technique to use residual estimates as refinement indicators [35]. For a critical point (uh , vh ) of I from the finite-dimensional space Xh,D × Xh , we use an estimate of the residual I ′ (uh , vh ; ϕ, ψ) for ϕ ∈ H1D0 , ψ ∈ H1 (Ω)
(3.1) ∗ in the H1D0 (Ω) × H1 (Ω) norm as a refinement indicator. Note that (3.1) is well defined since uh , vh ∈ W1,∞ (Ω). We derive an a posteriori upper bound on this residual in the first part of this section. In the second part we propose two adaptive finite element algorithms based on this bound to determine mesh refinement. 3.1. Residual Estimates. The following interpolation results will be needed for the subsequent residual estimate. Henceforth we use . to denote ≤ C where the constant C depends only on the shape regularity parameter K of the mesh but not on the mesh size. For a node with position xi , a face e, or an element T , let ωxi , ωe and ωT denote the union of all elements which touch the respective sets (or point). Further, let ∆i be the maximal ball contained in ωxi , if xi ∈ Ω, and the maximal ball intersected by Ω if i ∈ ∂Ω (cf. [36] for more detail). For χ ∈ H1 (Ω), define the quasi-interpolants Πh χ ∈ Xh and Πh,0 χ ∈ Xh,D0 as follows: Z Z X X χ dx λi (x). |∆i |−1 Πh χ (x) := χ dx λi (x), Πh,0 χ (x) := |∆i |−1 ∆i
i∈Nh
∆i
i∈Nh \Nh,D
(3.2) This quasi-interpolant was defined by Verf¨ urth [36] and was shown to satisfy the following approximation results. Let χ ∈ H1 (Ω) and Πh χ be as above; then, for all T ∈ Th and e ∈ Eh , kΠh χ − χkHs (T ) . h1−s T k∇χkL2 (ωT ) , kΠh χ − χkL2 (e) .
0 ≤ s ≤ 1,
h1/2 e k∇χkL2 (ωe ) .
(3.3) (3.4)
The same approximation result holds taking χ ∈ H1D0 (Ω) and replacing Πχ with Πh,0 χ. We also have the following approximation result for the nodal interpolant; see [13, Section 4.4]. For all T ∈ Th and w ∈ Ws,∞ (T ), kw − Ph wkL∞ (T ) . hsT |w|Ws,∞ (T ) ,
1 ≤ s ≤ 2.
(3.5)
Before stating the main proposition of this section it will be useful to introduce the following definition. For wh ∈ Xh and all e ∈ Eh , ∇wh |T1 − ∇wh |T2 , if e ∈ EI,h , with e = T¯1 ∩ T¯2 for some T1 , T2 ∈ Th , J∇wh Ke := |∇wh · n|e , if e ∈ EN,h ∪ ED,h , where n is the outer unit normal vector to ∂Ω. Proposition 3.1. Let uh ∈ Xh,D , vh ∈ Xh be such that Ih′ (uh , vh ; ϕh , ψh ) = 0 for all ϕh ∈ Xh,D0 , ψh ∈ Xh ; then |I ′ (uh , vh ; ϕ, ψ)| . µh k∇ϕkL2 (Ω) + νh k∇ψkL2 (Ω)
∀ϕ ∈ H1D0 (Ω), ∀ψ ∈ H1 (Ω),
where µh , νh are defined as follows: µh :=
h X
T ∈Th
|µT (uh , vh )|2
i1/2
,
νh :=
h X
T ∈Th
9
|νT (uh , vh )|2
i1/2
,
(3.6)
where, |µT (uh , vh )|2 := h4T k∇vh k4L∞ (T ) k∇uh k2L2 (T ) + h2T k2vh (∇vh · ∇uh )k2L2 (T ) X he kvh2 + ηk2L2 (e) J∇uh K2e , + e⊆∂T ∩(EI,h ∪EN,h )
2
|νT (uh , vh )| :=
h4T k∇vh k2L∞ (T ) k|∇uh |2 + ε2
X
e⊆∂T
+ αk2L2 (T ) + h2T k(|∇uh |2 + α)vh − αk2L2 (T )
he kJ∇vh Ke k2L2 (e) .
Proof. Since Ih′ (uh , vh ; ϕh , ψh ) = 0 for all ϕh ∈ Xh,D0 , ψ ∈ Xh it follows that ∀ϕh ∈ XhD0 , and
ah (vh ; uh , ϕh ) = 0
(3.7)
∀ψh ∈ Xh .
bh (uh ; vh , ψh ) = 0
(3.8)
Let us fix ϕ ∈ H1D0 (Ω) and ψ ∈ H1 (Ω); then |I ′ (uh , vh ; ϕ, ψ)| ≤ 2 |a(vh ; uh , ϕ)| + 2 |b(uh ; vh , ψ)|. Hence we shall examine the two functionals ϕ 7→ a(vh ; uh , ϕ) and ψ 7→ b(uh ; vh , ψ) separately. We begin by considering the former. By (3.7) we have |a(vh ; uh , ϕ)| ≤ |a(vh ; uh , ϕ − ϕh )| + |a(vh ; uh , ϕh ) − ah (vh ; uh , ϕh )|
∀ϕh ∈ Xh,D0 .
(3.9)
Examining the first term in (3.9), X Z |a(vh ; uh , ϕ − ϕh )| ≤ (vh2 + η)∇uh · ∇(ϕ − ϕh ) dx T ∈Th
T
X nZ = −2vh (∇vh · ∇uh )(ϕ − ϕh ) dx T
T ∈Th
+
Z
∂T
.
X
T ∈Th
o (vh2 + η)∇uh · n (ϕ − ϕh ) ds
k2vh (∇vh · ∇uh )kL2 (T ) kϕ − ϕh kL2 (T ) X
+
e⊆(EI,h ∪EN,h )
≤
X
T ∈Th
+
Z
e
J∇uh Ke |vh2 + η||ϕ − ϕh | ds
k2vh (∇vh · ∇uh )kL2 (T ) kϕ − ϕh kL2 (T ) X
e⊆(EI,h ∪EN,h )
J∇uh Ke kvh2 + ηkL2 (e) kϕ − ϕh kL2 (e) .
Since this inequality is true for all ϕh ∈ Xh,D0 , we choose ϕh = Πh,0 ϕ and use the approximation 10
results (3.3) and (3.4) to obtain |a(vh ; uh , ϕ − ϕh )| .
X
T ∈Th
hT k2vh (∇vh · ∇uh )kL2 (T ) k∇ϕkL2 (ωT ) X
+
e⊆(EI,h ∪EN,h )
.
h X
T ∈Th
+ .
h
2 h1/2 e kvh + ηkL2 (e) J∇uh Ke k∇ϕkL2 (ωe )
h2T k2vh (∇vh · ∇uh )k2L2 (T ) X
e⊆(EI,h ∪EN,h )
i1/2
he kvh2 + ηk2L2 (e) J∇uh K2e
h X
h2T k2vh (∇vh · ∇uh )k2L2 (T )
+
X
T ∈Th
k∇ϕkL2 (Ω)
e⊆∂T ∩(EI,h ∪EN,h )
i1/2
he kvh2 + ηk2L2 (e) J∇uh K2e
k∇ϕkL2 (Ω)
i1/2
k∇ϕkL2 (Ω) .
Using the approximation result (3.5) for the nodal interpolant and (3.3) with s = 1, we can bound the second term in (3.9) as follows: Z 2 2 |a(vh ; uh , ϕh ) − ah (vh ; uh , ϕh )| = {vh − Ph (vh )}∇uh · ∇ϕh dx Ω X Z ≤ |vh2 − Ph (vh2 )| |∇uh | |∇ϕh | dx T
T ∈Th
≤ .
X
kvh2 − Ph (vh2 )kL∞ (T ) k∇uh kL2 (T ) k∇ϕh kL2 (T )
X
h2T k∇vh k2L∞ (T ) k∇uh kL2 (T ) k∇ϕkL2 (ωT )
T ∈Th
T ∈Th
.
(3.10)
h X
T ∈Th
h4T k∇vh k4L∞ (T ) k∇uh k2L2 (T )
i1/2
k∇ϕkL2 (Ω) .
Therefore,
|a(vh ; uh , ϕ)| .
h X
T ∈Th
|µT (uh , vh )|2
i1/2
k∇ϕkL2 (Ω) ,
where |µT (uh , vh )|2 = h4T k∇vh k4L∞ (T ) k∇uh k2L2 (T ) + h2T k2vh (∇vh · ∇uh )k2L2 (T ) X he kvh2 + ηk2L2 (e) J∇uh K2e . + e⊆∂T ∩(EI,h ∪EN,h )
Now let us consider ψ 7→ b(uh ; vh , ψ). For all ψh ∈ Xh it follows from (3.8) that |b(uh ; vh , ψ)| ≤ |b(uh ; vh , ψ − ψh )| + |b(uh ; vh , ψh ) − bh (uh ; vh , ψh )|. 11
(3.11)
Considering the first term in (3.11) X Z |b(uh ; vh , ψ − ψh )| ≤ { (|∇uh |2 + α)vh − α}(ψ − ψh ) + ε∇vh · ∇(ψ − ψh ) dx T
T ∈Th
≤
X
T ∈Th
k(|∇uh |2 + α)vh − αkL2 (T ) kψ − ψh kL2 (T )
X Z + ε
∂T
T ∈Th
.
(ψ − ψh ) ∇vh · n ds
X
k ( |∇uh |2 + α)vh − αkL2 (T ) kψ − ψh kL2 (T )
+ε
X
T ∈Th
e⊆Eh
kψ − ψh kL2 (e) kJ∇vh Ke kL2 (e) .
This is true for all ψh ∈ Xh , so taking ψh = Πh ψ it follows that |b(uh ; vh , ψ − ψh )| .
X
T ∈Th
+ε
hT k ( |∇uh |2 + α)vh − αkL2 (T ) k∇ψkL2 (ωT )
X
h1/2 e kJ∇vh Ke kL2 (e) k∇ψkL2 (ωe )
e⊆Eh
.
h X
T ∈Th
+ε
h2T k( |∇uh |2 + α)vh − αk2L2 (T )
h X
e⊆Eh
.
h X T ∈Th
he kJ∇vh Ke k2L2 (e)
k∇ψkL2 (Ω)
k∇ψkL2 (Ω)
h2T k (|∇uh |2 + α)vh − αk2L2 (T )
X
+ ε2
i1/2
i1/2
e⊆∂T
he kJ∇vh Ke k2L2 (e)
i1/2
k∇ψkL2 (Ω) .
Finally we bound the second term in (3.11), noting that we have now fixed ψh = Πh ψ: Z 2 |b(uh ; vh , ψh ) − bh (uh ; vh , ψh )| = (vh ψh − Ph (vh ψh ))(|∇uh | + α) dx Ω X Z ≤ |vh ψh − Ph (vh ψh )| (|∇uh |2 + α) dx T
T ∈Th
.
N/2
X
hT
X
hT
X
hT
T ∈Th
.
k |∇uh |2 + αkL2 (T ) kvh ψh − Ph (vh ψh )kL∞ (T )
(N/2)+2
k |∇uh |2 + αkL2 (T ) |vh ψh |W2,∞ (T )
(N/2)+2
k |∇uh |2 + αkL2 (T ) k∇vh kL∞ (Ω) k∇ψh kL∞ (T ) .
T ∈Th
.
(3.12)
T ∈Th
Using the equivalence of norms on a finite-dimensional space we have |b(uh ; vh , ψh ) − bh (uh ; vh , ψh )| . ≤
X
(N/2)+2−(N/2)
hT
T ∈Th
h X
T ∈Th
k |∇uh |2 + αkL2 (T ) k∇vh kL∞ (Ω) k∇ψh kL2 (T )
h4T k|∇uh |2 + αk2L2 (T ) k∇vh k2L∞ (Ω) 12
i1/2
k∇ψkL2 (Ω) .
Therefore, |b(uh ; vh , ψ)| .
h X
T ∈Th
|νT (uh , vh )|2
i1/2
k∇ψkL2 (Ω) ,
where |νT (uh , vh )|2 = h4T k|∇uh |2 + αk2L2 (T ) |∇vh |2 + h2T k(|∇uh |2 + α)vh − αk2L2 (T ) X + ε2 he kJ∇vh Ke k2L2 (e) . e⊆∂T
We use µT (uh , vh ) as a local refinement indicator for uh , νT (uh , vh ) as a local refinement indicator for vh and define h 2 2 i1/2 ET (uh , vh ) := µT (uh , vh ) + νT (uh , vh ) for all T ∈ Th .
(3.13)
∗ Remark 3. We are in fact estimating the H1D0 (Ω) × H1 (Ω) norm of the gradient of I, since h X i1/2 |I ′ (uh , vh ; ϕ, ψ)| 2 . . |E (u , v )| kI ′ (uh , vh )k(H1D (Ω)×H1 (Ω))∗ = sup h T h h i1/2 0 ϕ∈H1 (Ω), 2 2 D0 T ∈T h + kψk kϕk H1 (Ω) H1 (Ω) 1 ψ∈H (Ω)
3.2. Adaptive Algorithm. We now propose two adaptive algorithms for computing local minimizers of I. The difference between the two algorithms is the stage at which the adaptive refinement of the mesh takes place. In Algorithm 1 this occurs after the alternate minimization algorithm has converged, whilst in Algorithm 2 the mesh is refined at each step of the alternate minimization algorithm. There are two user-specified tolerances associated with the algorithms: VTOL is the tolerance which determines when to halt the alternate minimization loop, whilst REFTOL is the tolerance determining when to halt the refinement loop. The marking parameter θ is a fixed number lying in the interval (0, 1]. In Algorithm 1 we denote the mesh at each level of refinement by Tj with hj = maxT ∈Tj diam(T ) for j ∈ N. In Algorithm 2 the mesh is also dependent on the alternate minimization step; accordingly, within each alternate minimization step n ∈ {m/2 : m ∈ N} we denote the mesh at each level of refinement by Tjn with hnj = maxT ∈Tjn diam(T ) for j ∈ N. Algorithm 1 1. Input: Initial crack field v0 and initial mesh T1 . 2. For j = 1, 2, . . . (a) Set vj1 = vj−1 ; (b) For n = 1, 2, . . . • unj := argmin {Ihj (z, vjn ) : z ∈ Xhj ,D }; • vjn+1 := argmin {Ihj (unj , z) : z ∈ Xhj }; • Repeat until kvjn+1 − vjn kL∞ (Ω) < VTOL.
(c) Set vj = vjn+1 , uj = unj ; P 2 1/2 > REFTOL, (d) If T ∈Tj |ET (uj , vj )|
P P • Find the smallest set Mj ⊆ Tj such that T ∈Mj |ET (uj , vj )|2 ≥ θ T ∈Tj |ET (uj , vj )|2 ; • Refine elements in Mj to obtain new mesh Tj+1 . 13
(e) Repeat until
P
T ∈Tj
| ET (uj , vj ) |2
3. Set uh (tk ) = uj and vh (tk ) = vj .
1/2
≤ REFTOL.
Algorithm 2 1. Input: Initial crack field v 1 and initial mesh T 1/2 . 2. For n = 1, 2, . . . (a) Set T1n = T n−1/2 (b) For j = 1, 2, . . . • Compute unj := argmin {Ihnj (z, v n ) : z ∈ Xhnj ,D }; √ P n n 2 1/2 > REFTOL/ 2, • If T ∈Tjn | µT (uj , v ) | P P • Find the smallest set Mj ⊆ Tjn such that T ∈Mj |µT (unj , vj )|2 ≥ θ T ∈T n |µT (unj , vj )|2 ; j n • Refine elements in Mj to obtain new mesh Tj+1 . √ P n n 2 1/2 ≤ REFTOL/ 2. • Repeat until T ∈T n | µT (uj , v ) | j
n
unj ,
n
(c) Set u = T = (d) For j = 1, 2, . . .
n+1/2
Tjn
and T1
= T n.
• Compute vjn+1 := argmin {Ihn+1/2 (un , z) : z ∈ Xhn+1/2 }; j j √ P n n+1 2 1/2 2, ) | > REFTOL/ • If n+1/2 | νT (u , v j T ∈T j
n+1/2
such that • Find the smallest set Mj ⊆ Tj P P n n+1 2 )| ≥ θ T ∈T n+1/2 |νT (un , vjn+1 )|2 ; T ∈Mj |νT (u , vj j
n+1/2
• Refine elements in Mj to obtain new mesh Tj+1 . √ P n n+1 2 1/2 )| ≤ REFTOL/ 2. • Repeat until n+1/2 | νT (u , v j T ∈T j
(e) Set v
n+1
=
vjn+1 and n+1
n+1/2
T n+1/2 = Tj
;
n
(f) Repeat until kv − v kL∞ (Ω) < VTOL. 3. Set vh (tk ) = v n+1 and uh (tk ) = un . The marking strategy used to identify elements for refinement is due to D¨ orfler [20]. In our implementation, the refinement of the mesh is achieved using the newest node bisection method [28], which guarantees a uniform bound on the shape-regularity of the generated meshes. Remark 4. We return to address a possible criticism, briefly mentioned at the end of Section 1.2, that a locally refined mesh may favour certain crack paths over other. For example, for Algorithm 1, it is quite conceivable that an ill-chosen initial mesh may have this effect. Note, however, that Steps (a) and (c) in Algorithm 2 in fact compute the exact solution of the alternating minimization steps, up to a prescribed tolerance. For example, suppose that the exact solution for one alternating minimization step initiates a new crack in a region of the domain where the mesh is coarse; then, the adaptive mesh refinement algorithm of Step (a) or (c) will be forced to refine the mesh in that region, provided the refinement tolerance is set sufficiently small. Thus, we strongly believe that Algorithm 2 is designed so that local adaptive refinement will not influence the formation of cracks beyond the usual numerical perturbation. In fact, local mesh refinement allows us to use a much smaller regularization parameters ǫ, and thus allows a more accurate computation of crack paths. Finally, we present a modification of Algorithm 2, which is useful for theoretical purposes. Algorithm 2′ : In step n of Algorithm 2′ , we replace REFTOL by REFTOLn , and we require that REFTOLn → 0 as n → ∞. Furthermore, we remove the termination condition (e). Remark 5. In order for Algorithms 2 and 2′ to be meaningful, steps (b) and (d) in Algorithm 2 need to terminate after a finite number of iterations. Without the mass lumping approximation, 14
this would follow immediately from standard convergence results for adaptive finite element methods for linear elliptic problems [14]. It is beyond the scope of this paper to prove that the same results remain true in the presence of mass lumping, however, we do not expect severe difficulties in extending existing results. Instead, in Theorem 4.2 we shall make the following assumption: Steps (a) and (c) in Algorithm 2 terminate in a finite number of iterations.
(A)
We now state some properties of the sequences generated by the preceding algorithms, which will prove useful later in the convergence analysis. Lemma 3.2. The sequence ((uj , vj ))∞ j=1 ⊆ Xhj ,D × Xhj generated by Algorithm 1 satisfies the following properties: 1. 0 ≤ vj (x) ≤ 1 on Ω for all j ∈ N; 1 1 2. ((uj , vj ))∞ j=1 is bounded in H (Ω) × H (Ω). Proof. The first property follows directly from the proof of Proposition 3. In order to prove the second property, note that uj := unj , vj := vjn+1 where ahj (vjn ; unj , ϕj ) = 0 bhj (unj ; vjn+1 , ψj )
=0
∀ϕj ∈ Xhj ,D0 ,
(3.14)
∀ψj ∈ Xhj .
(3.15)
Taking ϕj = unj − Phj g ∈ Xhj ,D0 in (3.14) we have Z (Phj ((vjn )2 ) + η)∇unj · ∇(unj − Phj g) dx = 0. Ω
Hence, ηk∇unj k2L2 (Ω) ≤
Z
Ω
(Phj ((vjn )2 ) + η)∇unj · ∇Phj g dx
≤ (1 + η)k∇unj kL2 (Ω) k∇Phj gkL2 (Ω) , and since we have assumed g ∈ Xhj ,D , k∇unj kL2 (Ω) ≤
1+η k∇gkL2 (Ω) . η
Therefore (k∇uj kL2 (Ω) )∞ j=1 is a bounded sequence, which by a variant of the Friedrichs inequality 1 implies that (kuj kL2 (Ω) )∞ j=1 is bounded as well; to see this, note that uj − Phj g ∈ HD0 (Ω), and so kuj kL2 (Ω) ≤ kuj − gkL2 (Ω) + kgkL2 (Ω) ≤ c∗ k∇(uj − g)kL2 (Ω) + kgkL2 (Ω) ≤ c∗ k∇uj kL2 (Ω) + (c2∗ + 1)1/2 kgkH1 (Ω) ,
where c∗ = c∗ (Ω) > 0 is the constant in the Friedrichs inequality. This implies that (uj )∞ j=1 is bounded in H1 (Ω). Taking ψj = vjn+1 ∈ Xj in (3.15) we have Z Z Z |∇vjn+1 |2 dx = 0, vjn+1 (vjn+1 − 1) dx + ε Phj ((vjn+1 )2 )|∇unj |2 dx + α Ω
Ω
Ω
which implies that εk∇vjn+1 k2L2 (Ω) ≤ α
Z
Ω
(1 − vjn+1 )vjn+1 dx ≤
α|Ω| . 4
2 ∞ ∞ Hence (k∇vj k)∞ j=1 is bounded in L (Ω), and since (vj )j=1 is bounded in L (Ω) it is also bounded 2 ∞ 1 1 in L (Ω). Thus ((uj , vj ))j=1 is bounded in H (Ω) × H (Ω).
15
′ Lemma 3.3. The sequence (un , v n )∞ n=1 ⊆ Xhn ,D × Xhn generated by Algorithm 2 satisfies the following properties: 1. 0 ≤ v n (x) ≤ 1 on Ω for all n ∈ N; 1 1 2. ((un , v n ))∞ n=1 is a bounded sequence in H (Ω) × H (Ω); n n n n n n n 3. |I(u , v ) − Ih (u , v )| . νh k∇v kL2 (Ω) for all n ∈ N; 4. lim inf n→∞ I(un , v n ) ≤ lim inf n→∞ I(un−1 , v n ) ≤ lim inf n→∞ I(un−1 , v n−1 ). Proof. Properties 1 and 2 follow in a similar manner to that of Lemma 3.2. In order to show Property 3 note that Z |I(un , v n ) − Ihn (un , v n )| = [(v n )2 − Phn ((v n )2 )](|∇un |2 + α) dx Ω
= |b(un ; v n , v n ) − bhn (un ; v n , v n )|. Thus using the derivation of the residual estimate it follows that |I(un , v n ) − Ihn (un , v n )| . νhn k∇v n kL2 (Ω) . To show the left-hand inequality in Property 4 let u ˜ := argminXhn I(·, v n ). Then, I(un , v n ) ≤ Ihn (un , v n ) + |I(un , v n ) − Ihn (un , v n )|
u, v n ) + |I(un , v n ) − Ihn (un , v n )| ≤ Ihn (˜ u, v n )| + |I(un , v n ) − Ihn (un , v n )| ≤ I(˜ u, v n ) + |I(˜ u, v n ) − Ihn (˜
u, v n )| + |I(un , v n ) − Ihn (un , v n )|. ≤ I(un−1 , v n ) + |I(˜ u, v n ) − Ihn (˜
Since limn→∞ REFTOLn = 0 it follows that limn→∞ νhn+1/2 = 0. Thus using Property 3 and letting n → ∞ we have lim inf I(un , v n ) ≤ lim inf I(un−1 , v n ). n→∞
n→∞
The right-hand inequality follows in a similar manner. 4. Convergence Analysis. In this section we state and prove two results which support use of the adaptive algorithms proposed in the previous section. Provided that Algorithm 1 terminates for any given input (we are in fact unable to prove this at present) then Theorem 4.1 demonstrates convergence of the numerical solutions for decreasing tolerance REFTOL to a critical point of I. Theorem 4.1. Assume that Ω is an open bounded domain in RN . Suppose that there exists 1 1 ∞ a sequence ((uj , vj ))∞ j=1 ⊂ HD (Ω) × (H (Ω) ∩ L (Ω)), 0 ≤ vj (x) ≤ 1 a.e. in Ω, such that a(vj ; uj , ϕ) ≤ µj k∇ϕkL2 (Ω) b(uj ; vj , ψ) ≤ νj kψkH1 (Ω)
∀ϕ ∈ H1D0 (Ω), 1
(4.1) ∞
∀ψ ∈ H (Ω) ∩ L (Ω),
(4.2)
for some µj , νj ∈ R≥0 with µj , νj → 0 as j → ∞. Suppose also that ((uj , vj ))∞ j=1 is a bounded 1 1 sequence in H (Ω) × H (Ω). 1 1 Then, there exists a subsequence of ((uj , vj ))∞ j=1 (not relabelled) and (u, v) ∈ HD (Ω) × H (Ω), 1 1 with 0 ≤ v(x) ≤ 1 a.e. in Ω, such that uj → u in H (Ω) and vj → v in H (Ω) as j → ∞. Moreover, u and v satisfy a(v; u, ϕ) = 0 b(u; v, ψ) = 0
∀ϕ ∈ H1D0 (Ω), 1
(4.3) ∞
∀ψ ∈ H (Ω) ∩ L (Ω),
(4.4)
that is, (u, v) is a critical point of I in H1D (Ω) × H1 (Ω) ∩ L∞ (Ω). Theorem 4.2 proves that the sequence (un , v n ) computed by the Algorithm 2′ (which is designed without termination criterion) converges to a critical point of I and thus, subject to a justification of Assumption A, establishes 16
the convergence of the Algorithm 2′ . To the best of our knowledge this is the first convergence result for an adaptive finite element algorithm for a nonconvex minimization problem. Theorem 4.2. Assume that Ω is an open bounded domain in RN . Let ((un , vn ))∞ n=1 ⊂ 1 HD (Ω) × H1 (Ω) be the sequence generated by Algorithm 2′ under assumption (A). ∞ 1 1 Then, there exists a subsequence ((unj , vnj ))∞ j=1 of ((un , vn ))n=1 and (u, v) ∈ HD (Ω) × H (Ω), 1 1 with 0 ≤ v(x) ≤ 1 a.e. in Ω, such that unj → u in H (Ω) and vnj → v in H (Ω) as j → ∞. In addition, u and v satisfy a(v; u, ϕ) = 0 b(u; v, ψ) = 0
∀ϕ ∈ H1D0 (Ω),
(4.5)
1
∞
∀ψ ∈ H (Ω) ∩ L (Ω),
(4.6)
that is, (u, v) is a critical point of I in H1D (Ω) × H1 (Ω) ∩ L∞ (Ω). Proof. [of Theorem 4.1] Step 1. Existence of a convergent subsequence of ((uj , vj ))∞ j=1 : 1 1 The sequence ((uj , vj ))∞ is bounded in H (Ω) × H (Ω). Since H1 (Ω) is a Hilbert space, there j=1 ∞ exists a subsequence of ((uj , vj ))j=1 (not relabelled) such that (uj , vj ) ⇀ (u, v)
in H1 (Ω) × H1 (Ω) as j → ∞,
As H1D (Ω) is a convex and closed subset of H1 (Ω) it is also weakly closed (c.f. Proposition 2.5 on p.35 of [6]). Hence, we have that u ∈ H1D (Ω). Similarly, as the set K := {w ∈ H1 (Ω) : 0 ≤ w(x) ≤ 1 a.e. x ∈ Ω} is a convex closed subset of H1 (Ω), and since vj ∈ K for all j ∈ N, it follows that 0 ≤ v(x) ≤ 1 a.e. in Ω. R Step 2. limj→∞ Ω |v − vj | |∇w|2 dx = 0 for all w ∈ H1 (Ω): (Recall that we have not relabelled the convergent subsequence.) This result will prove repeat1 ∞ edly useful in the subsequent R analysis. Fix2 w ∈ H (Ω) and let R (vjk )k=1 be a2 subsequence of ∞ (vj )j=1 such that limk→∞ Ω |vjk − v| |∇w| dx = lim supj→∞ Ω |vj − v| |∇w| dx, and vjk → v a.e. in Dominated ConvergenceR Theorem [37, Section 5.2], we have R Ω. Using Lebesgue’s 2 |∇w| dx = 0 and hence lim supj→∞ Ω |vj − v| |∇w|2 dx = 0. It thus follimk→∞ Ω |vjk − v| R 2 lows that limj→∞ Ω |vj − v| |∇w| dx = 0.
Step 3. a(v; u, ϕ) = 0 for all ϕ ∈ H1D0 (Ω): Fixing ϕ ∈ H1D0 (Ω) we have a(v; u, ϕ) = =
Z
ZΩ
(v 2 + η)∇u · ∇ϕ dx
(v 2 + η)∇(u − uj ) · ∇ϕ dx + Ω Z + (v 2 − vj2 )∇uj · ∇ϕ dx
Z
Ω
(vj2 + η)∇uj · ∇ϕ dx
Ω
=: Rj + Sj + Tj . Our aim is to show that Rj , Sj , Tj → 0 as j → ∞. First, consider Rj =
Z
Ω
(v 2 + η)∇ϕ · ∇(u − uj ) dx.
2 N We know that (v 2 +η)∇ϕ ∈ (L2 (Ω))N ; so, by the weak convergence of (∇uj )∞ j=1 to ∇u in (L (Ω)) , we obtain that Rj → 0 as j → ∞. Turning our attention to Sj , (4.1) implies that
|Sj | ≤ µj k∇ϕkL2 (Ω) → 0 17
as j → ∞.
Finally, we estimate Tj by Z
|v + vj ||v − vj ||∇uj ||∇ϕ| dx Z 1/2 k∇uj kL2 (Ω) , ≤2 |v − vj | |∇ϕ|2 dx
|Tj | ≤
Ω
Ω
2 N where we used the fact |v + vj | ≤ 2 and |v − vj | ≤ 1. Since (∇uj )∞ j=1 is bounded in (L (Ω)) , and using Step 2 with w = ϕ, we have Tj → 0 as j → ∞. Thus we conclude that a(v; u, ϕ) = 0 for all ϕ ∈ H1D0 (Ω).
Step 4. ∇uj → ∇u strongly in (L2 (Ω))N : Note that Z 2 ηk∇u − ∇uj kL2 (Ω) ≤ (vj2 + η)∇(u − uj ) · ∇(u − uj ) dx Ω Z Z = − (vj2 + η)∇uj · ∇(u − uj ) dx + (vj2 + η)∇u · ∇(u − uj ) dx Ω Ω Z 2 ≤ µj k∇(u − uj )kL2 (Ω) + (vj + η)∇u · ∇(u − uj ) dx. Ω
Since u − uj ∈
H1D0 (Ω),
we have a(v; u, u − uj ) = 0, and it follows that Z ηk∇u − ∇uj k2L2 (Ω) ≤ µj k∇(u − uj )kL2 (Ω) + (vj2 − v 2 )∇u · ∇(u − uj ) dx, Ω
≤ µj k∇(u − uj )kL2 (Ω) + k(vj2 − v 2 )∇ukL2 (Ω) k∇(u − uj )kL2 (Ω) .
Dividing by k∇u − ∇uj k, and using the fact that |vj2 − v 2 |2 ≤ 4|vj − v|2 ≤ 4|vj − v|, we obtain Z 1/2 2 2 ηk∇u − ∇uj kL2 (Ω) ≤ µj + k(vj − v )∇ukL2 (Ω) ≤ µj + 2 |vj − v| |∇u|2 dx . Ω
Since µj → 0 as j → ∞, and using Step 2, we deduce that limj→∞ k∇u − ∇uj kL2 (Ω) = 0. Step 5. b(u; v, ψ) = 0 for all ψ ∈ H1 (Ω) ∩ L∞ (Ω): Fixing ψ ∈ H1 (Ω) ∩ L∞ (Ω), we have Z (|∇u|2 + α)vψ − αψ + ε∇v · ∇ψ dx b(u; v, ψ) = ZΩ (|∇u|2 + α)(v − vj )ψ + ε∇(v − vj ) · ∇ψ dx = Ω Z (|∇uj |2 + α)vj ψ − αψ + ε∇vj · ∇ψ dx + ZΩ (|∇u|2 − |∇uj |2 )vj ψ dx + Ω
=: Rj + Sj + Tj .
We can estimate Rj by Z |Rj | ≤ |∇u|2 (v − vj )ψ dx + αkv − vj kL2 (Ω) kψkL2 (Ω) Ω Z + ε ∇(v − vj ) · ∇ψ dx Ω Z |∇u|2 |v − vj | dx + αkv − vj kL2 (Ω) kψkL2 (Ω) ≤ kψkL∞ (Ω) Ω Z +ε ∇(v − vj ) · ∇ψ dx . Ω
18
Since vj → v in L2 (Ω), ∇vj ⇀ ∇v in (L2 (Ω))N , and using Step 2 with w = u, we deduce that Rj → 0 as j → ∞. Now, from (4.2) we know that |Sj | ≤ ηj kψkH1 (Ω) → 0, and therefore Sj → 0 as j → ∞. Lastly, we estimate Tj by Z |Tj | = (|∇u|2 − |∇uj |2 )vj ψ dx Ω Z ≤ kψkL∞ (Ω) |∇(u − uj )||∇(u + uj )| dx Ω
≤ kψkL∞ (Ω) k∇(u − uj )kL2 (Ω) k∇(u + uj )kL2 (Ω) . Since ∇uj → ∇u in (L2 (Ω))N , we obtain Tj → 0 as j → ∞. Hence b(u; v, ψ) = 0 for all ψ ∈ H1 (Ω) ∩ L∞ (Ω). Step 6. ∇vj → ∇v strongly in (L2 (Ω))N : Observe that Z (|∇u|2 + α)|v − vj |2 + ε(∇v − ∇vj ) · (∇v − ∇vj ) dx εk∇v − ∇vj k2L2 (Ω) ≤ Ω Z (|∇u|2 + α)vj (v − vj ) + ε∇vj · ∇(v − vj ) − α(v − vj ) dx =− ZΩ (|∇u|2 + α)v(v − vj ) + ε∇v · ∇(v − vj ) − α(v − vj ) dx + Ω Z = −b(uj ; vj , v − vj ) + (|∇uj |2 − |∇u|2 )vj (v − vj ) dx + b(u; v, v − vj ). Ω
We have already shown that b(u; v, v − vj ) = 0, and hence, Z 2 |∇uj + ∇u||∇uj − ∇u||vj ||v − vj | dx εk∇v − ∇vj kL2 (Ω) ≤ |b(uj ; vj , v − vj )| + Ω
≤ νj kv − vj kH1 (Ω) + k∇u + ∇uj kL2 (Ω) k∇uj − ∇ukL2 (Ω) . ∞ Since νj → 0, ∇uj → ∇u in (L2 (Ω))N , and since (k∇u + ∇uj kL2 (Ω) )∞ j=1 and (kv − vj kH1 (Ω) )j=1 are bounded, we conclude that ∇vj → ∇v in (L2 (Ω))N as j → ∞, which completes the proof of the theorem. Proof. [of Theorem 4.2] Step 1. Existence of convergent subsequences of (un , vn )∞ n=1 : 1 1 Since ((un , vn ))∞ is a bounded sequence in H (Ω) × H (Ω), it follows analogously to Step 1 in n=1 D ∞ the proof of Theorem 4.1 that there exists a subsequence (unj , vnj )∞ j=1 of (un , vn )n=1 such that
(unj , vnj ) ⇀ (u, v) in H1 (Ω) × H1 (Ω) as j → ∞, for some u ∈ H1D (Ω) and v ∈ H1 (Ω), with 0 ≤ v(x) ≤ 1 a.e. in Ω. It will be seen that careful labelling of the subsequence is vital in the proof. Upon extracting a further subsequence we may assume, without loss of generality, that vnj +1 ⇀ v ′ weakly in H1 for some v ′ ∈ H1 (Ω). Step 2. a(v; u, ϕ) = 0 for all ϕ ∈ H1D0 (Ω): By definition of (unj , vnj )∞ j=1 we have a(vnj ; unj , ϕ) ≤ µnj k∇ϕkL2 (Ω) for all j ∈ N and ϕ ∈ H1D0 (Ω), with µnj → 0 as j → ∞. It follows analogously to Steps 3 and 4 of the proof of Theorem 4.1 that a(v; u, ϕ) = 0 for all ϕ ∈ H1D0 (Ω), 19
and also that ∇unj → ∇u in (L2 (Ω))N as j → ∞. Step 3. Strong convergence of (vnj +1 )∞ j=1 : By the compact embedding of H1 (Ω) in L2 (Ω) it follows that vnj +1 → v ′ strongly in L2 (Ω). To obtain strong convergence in H1 (Ω), note that the sequence (unj , vnj +1 )∞ j=1 satisfies b(unj ; vnj +1 , ψ) ≤ νnj k∇ψkL2 (Ω) for all j ∈ N and for all ψ ∈ H1 (Ω). Taking ψ = (vnj+1 − v ′ ) ∈ H1 (Ω) it follows that Z Z vnj +1 (vnj +1 − v ′ ) |∇unj |2 dx + α (1 − vnj +1 )(vnj +1 − v ′ ) dx νnj k∇ψkL2 (Ω) ≥ Ω Ω Z +ε ∇vnj +1 · (∇vnj +1 − ∇v ′ ) dx, Ω
and therefore, Z Z Z ε |∇vnj +1 |2 dx ≤ vnj +1 (v ′ − vnj +1 )|∇u|2 dx + vnj +1 (v ′ − vnj +1 )(|∇unj |2 − |∇u|2 ) dx Ω Ω Ω Z Z ′ + α (1 − vnj +1 )(v − vnj +1 ) dx + ε ∇vnj +1 · ∇v ′ dx Ω
≤
Ω
+ νnj k∇ψkL2 (Ω) ,
Z
|v ′ − vnj +1 | |∇u|2 dx + k∇unj − ∇ukL2 (Ω) k∇unj + ∇ukL2 (Ω) Ω Z Z |v ′ − vnj +1 | dx + ε +α ∇vnj +1 · ∇v ′ dx + νnj k∇ψkL2 (Ω) . Ω
Ω
Letting j → ∞, we obtain
lim sup ε j→∞
Z
Ω
2
|∇vnj +1 | ≤ ε
Z
Ω
|∇v ′ |2 dx.
Using weak lower semi-continuity of the L2 -norm, we conclude that Z Z Z Z ′ 2 2 2 |∇v | ≤ lim inf |∇vnj +1 | dx ≤ lim sup |∇vnj +1 | dx ≤ |∇v ′ |2 dx. Ω
j→∞
j→∞
Ω
Ω
Ω
Thus we deduce that k∇vnj +1 kL2 (Ω) → k∇v ′ kL2 (Ω) , which, together with weak convergence, implies the strong convergence of vnj+1 to v ′ in H1 (Ω) as j → ∞. Step 4. v = v ′ ; b(u; v, ψ) = 0 for all ψ ∈ H1 (Ω): Since the sequence ((unj , vnj +1 )∞ j=1 satisfies b(unj ; vnj +1 , ψ) ≤ νnj k∇ψkL2 (Ω) for all j ∈ N and for all ψ ∈ H1 (Ω). Letting j → ∞, it follows analogously to Step 5 in the proof of Theorem 4.1, that b(u; v ′ , ψ) = 0 for all ψ ∈ H1 (Ω), and we are left to show that v = v ′ . Using the weak lower-semicontinuity of I together with Property 4 of Lemma 3.3 it follows that I(u, v) ≤ lim inf I(unj+1 , vnj+1 ) ≤ lim I(unj , vnj +1 ) = I(u, v ′ ). j→∞
j→∞
Since v ′ is a critical point of the strictly convex map v 7→ I(u, v), v ′ is its unique minimizer. Since I(u, v) = I(u, v ′ ) it follows that v = v ′ and hence (u, v) is a critical point of I. 20
5. A Computational Example. The purpose of this section is to present some numerical results for the adaptive algorithms proposed in Section 3.2. We will consider just one example chosen to address some of the relevant questions surrounding the method. We would like to stress the preliminary nature of the results and assert that there are still aspects of the results that we are yet to fully understand. It is our intention to explore these more fully in future work, together with an extension of the implementation to a more general setting. 5.1. A Two-Dimensional Anti-Plane Strain Example. The example chosen takes Ω to be the two-dimensional rectangular domain (0, 2) × (0, 2.2), containing a slit along {1} × [1.5, 2.2], and a circular hole with centre (0.3, 0.3), radius 0.2. This is shown in Figure 5.1 where the shaded part (0, 1) ∪ (1, 2) × (2, 2.2) is ΩD . The incremental anti-plane displacement g(x, t) is given as follows: −t, on (0, 1) × (2, 2.2), g(x, t) = t, on (1, 2) × (2, 2.2), where t = 0.01s for s = 0, 1, . . . , 150. In all subsequent experiments we fix η = 1 × 10−5 , VTOL = 1 × 10−3 , CRTOL = 1 × 10−4 and θ = 0.3. The initial triangulation of the domain, shown in Figure 5.1, is generated using the software package Triangle [31, 32]. At each time step the initial crack field v is taken to be the final computed v from the previous time step, with the exception of the first time step where it is taken to be v ≡ 1. In the implementation we do not enforce condition (2.2). 2.2 2.0
1.5
0.3 0 0 0.1 0.3 0.5
1.0
2.0
(a) Domain
(b) Initial mesh
Fig. 5.1. Computational domain
We first use Algorithm 2 to compute the evolution of the crack with ε = 0.02 and REFTOL = 0.05. Figure 5.2 shows the crack field variable v and the mesh at the final time step, together with the change in bulk, surface and total energies over time. The body first remains unfractured until time s ∼ 25 when the applied displacement becomes sufficiently large to initiate a crack at the slit tip. A period of slow continuous crack growth then follows. Finally at time s = 124 failure occurs, when the crack length rapidly increases at a single time splitting the body into two pieces. The sudden decrease in the total energy at this time reveals that the algorithm has computed local minimizers of the energy. In computing the evolution of the crack those elements of the initial mesh located near the crack path are repeatedly refined. The number of elements in the final mesh is 508947 with the smallest and largest elements having diameters of order 10−4 and 10−1 respectively. It is natural to consider how the choice of ε and REFTOL affect the evolution of the crack. Figure 5.3 shows the final crack path, computed for a range of parameter values using Algorithm 2. The values have been selected to demonstrate the wide variation of crack paths possible if the parameters are not chosen sufficiently small. In particular we believe these results highlight the importance of choosing ε suitably small, which in turn requires the use of a sufficiently fine mesh, close to the crack, to fully resolve the solution. This is clearly more easily achievable using an adaptively refined mesh than a uniform mesh. 21
2.5 Bulk Energy Surface Energy Total Energy 2
1.5
1
0.5
0 0
(a) Crack path: plot of v
(b) Final mesh
20
40
60
80
100
120
140
(c) Bulk and surface energy
Fig. 5.2. The evolution of the body computed by Algorithm 2 for ε = 0.02 and REFTOL = 0.05
(a) ε = 0.02, REFTOL = 0.3 (b) ε = 0.02, REFTOL = 0.2 (c) ε = 0.02, REFTOL = 0.1
(d) ε = 0.05, REFTOL = 0.5 (e) ε = 0.05, REFTOL = 0.1 (f) ε = 0.05, REFTOL = 0.05 Fig. 5.3. A comparison of the final crack path computed by Algorithm 2 for different values of ε and REFTOL
Thus far we have only considered Algorithm 2. We now compare the results of Algorithm 1 and Algorithm 2 with the choice of parameters ε = 0.02 and REFTOL = 0.05. A comparison by eye shows that the final crack path for Algorithm 1, shown in Figure 5.4, is very similar to that computed by Algorithm 2 (Figure 5.2 (a)). The failure time for Algorithm 1 is at s = 125, one time step later than that of Algorithm 2. In this experiment we observe that Algorithm 1 takes more alternate minimization steps than Algorithm 2; however, Algorithm 2 is required to take more refinement iterations and thus generates meshes with a greater number of elements. Consequently, while Algorithm 1 solves a greater number of linear systems than Algorithm 2, the linear systems solved in Algorithm 2 are generally larger. Finally, we note that in all the computations we observe that 0 ≤ vh (x) ≤ 1 for all x ∈ Ω, despite the fact we do not explicitly enforce condition (2.2) through regularity conditions on the mesh.
22
Fig. 5.4. The final crack path computed using Algorithm 1 for ε = 0.02 and REFTOL = 0.05
6. Conclusion. We have presented two adaptive algorithms for computing finite element approximations of local minimizers of the Ambrosio–Tortorelli functional. In this paper we have primarily focused on proving convergence results for the two algorithms. We have been able to show that, provided the associated residuals converge to zero, both algorithms generate a sequence of numerical solutions that converge to a critical point of I (as the termination tolerances go to zero). Our preliminary computational results demonstrate the potential of using this method. In particular, we believe that the algorithms enable us to accurately and reliably compute the evolution of the crack path using considerably fewer elements compared to simulations with uniform meshes. However, we have observed that the crack path can be sensitive to the choice of the parameters in the adaptive algorithm. We will investigate this further through more extensive testing of the algorithms in future work. We believe that the results presented are easily extendable to the cases of planar and threedimensional elasticity. We are currently working on extending the theory and implementation to these models. REFERENCES [1] L. Ambrosio, A compactness theorem for a new class of functions of bounded variation, Boll. Un. Mat. Ital. B (7), 3 (1989), pp. 857–881. [2] , Existence theory for a new class of variational problems, Arch. Rational Mech. Anal., 111 (1990), pp. 291–322. [3] L. Ambrosio and V. Tortorelli, Approximation of functionals depending on jumps by elliptic functionals via Γ-convergence, Comm. Pure. Appl. Math., 43 (1990), pp. 999–1036. [4] , On the approximation of free discontinuity problems, Boll. Un. Mat. Ital. B (7), 6 (1992), pp. 105–123. [5] N. Ambrosio, L. Fusco and D. Pallara, Functions of Bounded Variation and Free Discontinuity Problems, Oxford University Press, 2000. [6] V. Barbu and T. Precopeanu, Convexity and Optimization in Banach Spaces, Springer Verlag, 1986. [7] G. Bellittini and A. Coscia, Discrete approximation of a free discontinuity problem, Numer. Funct. Anal. Optim., 15 (1994), pp. 105–123. [8] B. Bourdin, Numerical implementation of the variational formulation for quasi-static brittle fracture, Interfaces Free Bound., 9 (2007), pp. 411–430. [9] , The variational formulation of brittle fracture: numerical implementation and extensions, IUTAM Symposium on Discretization Methods for Evolving Discontinuities, (2007). [10] B. Bourdin and A. Chambolle, Implementation of an adaptive finite element approximation of the Mumford-Shah functional, Numer. Math., 85 (2000), pp. 609–646. [11] B. Bourdin, G. Francfort, and J.-J. Marigo, Numerical experiments in revisited brittle fracture, J. Mech. Phys. Solids, 48 (2000), pp. 797–826. [12] A. Braides, Γ-Convergence for Beginners, vol. 22 of Oxford Lecture Series in mathematics and its Applications, Oxford University Press, 2002. [13] S. Brenner and L. Scott, The Mathematical Theory of Finite Element Methods, Springer Verlag, Second ed., 1994. [14] J. Cascon, C. Kreuzer, R. Nochetto, and K. Siebert, Quasi-optimal convergence rate for an adaptive finite element method. to appear in SIAM J. Numer. Anal. [15] A. Chambolle and G. Dal Maso, Discrete approximation of the Mumford-Shah functional in dimension 23
two, M2AN Math. Model. Numer. Anal., 33 (2000), pp. 609–646. [16] P. Ciarlet and P.-A. Raviart, Maximum principle and uniform convergence for the finite element method, Comput. Methods Appl. Mech. Engrg., 2 (1973), pp. 17–31. [17] G. Dal Maso, G. Francfort, and R. Toader, Quasistatic crack growth in nonlinear elasticity, Arch. Ration. Mech. Anal., 176 (2005). [18] G. Dal Maso and R. Toader, A model for the quasi-static growth of brittle fractures: existence and approximation results, Arch. Ration. Mech. Anal., 162 (2002), pp. 101–135. [19] M. De Giorgi, E. Carriero and A. Leaci, Existence theorem for a minimum problem with free discontinuity set, Rat. Mech. Anal., 108 (1989), pp. 195–218. ¨ rfler, A convergent adaptive algorithm for Poisson’s equation, SIAM J. Numer. Anal., 33 (1996), [20] W. Do pp. 1106–1124. [21] L. Evans and R. Gariepy, Measure Theory and Fine Properties of Functions, Studies in Advanced Mathematics, CRC Press, 1992. [22] M. Focardi, On the variational approach of free discontinuity problems in the vectorial case, Math. Models Methods Appl. Sci., 11 (2001), pp. 663–684. [23] G. Francfort and C. Larsen, Existence and convergence for quasi-static evolution in brittle fracture, Comm. Pure Appl. Math., 56 (2003), pp. 1465–1500. [24] G. Francfort and J.-J. Marigo, Revisiting brittle fracture as an energy minimization problem, J. Mech. Phys. Solids, 46 (1998), pp. 1319–1342. [25] A. Giacomini, Ambrosio-Tortorelli approximation of quasi-static evolution of brittle fractures, Calc. Var. Partial Differential Equations, 22 (2005), pp. 129–172. [26] A. Griffith, The phenomena of rupture and flow in solids, Philosophical Transactions of the Royal Society of London, 221 (1921), pp. 163–198. ˇ´ıˇ ¨ ki, Weakened acute type condition for tetrahedral triangulations [27] S. Korotov, M. Kr zek, and P. Neittaanma and the discrete maximum principle, Math. Comp., 70 (2001), pp. 107–119 (electronic). [28] W. Mitchell, A comparison of adaptive refinement techniques for elliptic problems, ACM Trans. Math. Softw., 15 (1989), pp. 326–247. [29] D. Mumford and J. Shah, Optimal approximation by piecewise smooth functions and associated variational problems, Comm. Pure Appl. Math., 42 (1989), pp. 577–685. [30] M. Negri and M. Paolini, Numerical minimization of the Mumford-Shah functional, Calcolo, 38 (2001), pp. 67–84. [31] J. Shewchuk, Triangle: Engineering a 2D quality mesh generator and delaunay triangulator, First Workshop on Applied Computational Geometry Proceedings, (1996), pp. 124–133. [32] , Delaunay refinement algorithms for triangular mesh generation, Computational Geometry: Theory and Applications, 22 (2002), pp. 21–74. [33] G. Strang and G. Fix, An analysis of the finite element method, Prentice-Hall Inc., 1973. [34] V. Thom´ ee, Galerkin finite element methods for parabolic problems, vol. 1054 of Lecture Notes in Mathematics, Springer-Verlag, 1984. ¨ rth, A Review of A Posteriori Error Estimation and Adaptive Mesh-Refinement Techniques, Wiley [35] R. Verfu & Teubner, 1996. ¨ rth, Error estimates for some quasi-interpolation operators, Mod´ [36] R. Verfu el. Math. Anal. Num´ er, 33 (1999), pp. 695–713. [37] A. Weir, Lebesgue integration and measure, Cambridge University Press, 1973.
24
