Approximation of fracture energies with $p$-growth via piecewise affine finite elements

The modeling of fracture problems within geometrically linear elasticity is often based on the space of generalized functions of bounded deformation $GSBD^p(\Omega)$, $p\in(1,\infty)$, their treatment is however hindered by the very low regularity of those functions and by the lack of appropriate density results. We construct here an approximation of $GSBD^p$ functions, for $p\in(1,\infty)$, with functions which are Lipschitz continuous away from a jump set which is a finite union of closed subsets of $C^1$ hypersurfaces. The strains of the approximating functions converge strongly in $L^p$ to the strain of the target, and the area of their jump sets converge to the area of the target. The key idea is to use piecewise affine functions on a suitable grid, which is obtained via the Freudhental partition of a cubic grid.


Introduction
The modeling of plasticity and fracture in a geometrically linear framework leads to vectorial variational problems in which the local energy depends on the symmetric part of the deformation gradient, and the deformation can jump in a set of finite (n − 1)-dimensional measure [22,28,29]. If one assumes that the total variation of the distributional symmetrized gradient is controlled by the energy then one deals with functions of bounded deformation, which are defined as the functions u ∈ L 1 (Ω; R n ) such that the distributional strain Eu := 1 2 (Du + Du T ) is a bounded measure [1,3,27,29,30].
Here Ω ⊂ R n is an open set, and BD(Ω) denotes the set of functions of bounded deformation on Ω.
In fracture problems one often deals with the proper subspace SBD p (Ω), which is characterized by the fact that the distributional strain Eu is the sum of an elastic part e(u)L n Ω, with e(u) ∈ L p (Ω; R n×n sym ), and a singular part [u] ν u H n−1 J u concentrated on a (n − 1)-rectifiable set of finite (n − 1)-dimensional measure, with ν u the approximate normal and [u] the jump of the traces of u across J u (see [4,5,8,22]). Typical fracture models, such as Griffith's model, do not, however, give control of the amplitude of the jump of u over the discontinuity set; a typical energy takes the form Ω f (e(u))dx + H n−1 (J u ), (1.1) which is the natural vectorial generalization of the Mumford-Shah functional in linear elasticity. The function f is assumed to be convex and to have p-growth at infinity. One is then lead to compactness results in the space GSBD p (Ω), which was introduced by Dal Maso in [20] and recalled in Section 2.
In the study of problems modeled in SBD p or GSBD p (see for example [25]) it is crucial to have good approximation results for functions in those spaces. On the one hand, smooth functions are dense in BD(Ω) in the strict topology, which entails weak convergence of the distributional strains. This is clearly not enough to ensure continuity of the energy in (1.1) along such approximating sequences. Indeed, the smooth approximants (which are typically obtained by mollification) replace both the discontinuities and the L p strain e(u) by smooth components, mixing fracture and elastic deformation. It is apparent that this will, in general, increase significantly the energy. In the scalar case from the point of view of applications to fracture mechanics, the functional setting for the problem is provided by SBV p (Ω; R n ) functions, and a density result which guarantees separate convergence of the two terms in (1.1) was obtained by Cortesani and Toader [19]. The approximants are still discontinuous, but the jump set has become regular (a finite union of simplexes) and each function is regular away from the jump set. The gradients converge strongly away from the jump set, and the jumps and the orientation of the jump set converge. More generally, in such a restricted framework one can allow the domain and the codomain to have different dimensions, in what follows we shall limit to comment the case of interest in this paper.
In the vector-valued case, and restricting to energies with quadratic growth, density of regular functions in SBD 2 (Ω) ∩ L 2 (Ω; R n ) was proven by Chambolle in 2004 for n = 2 [9] and then for n ≥ 3 [10]. His proof was extended to GSBD 2 (Ω) ∩ L 2 (Ω; R n ) by Iurlano [26]. Their result shows that any u ∈ GSBD 2 (Ω) ∩ L 2 (Ω; R n ) can be approximated by functions which are continuous away from a finite union of closed pieces of C 1 hypersurfaces, are Lipschitz continuous away from this set, with strong convergence of the strains and, in an appropriate sense, of the discontinuities. This permits to obtain convergence of energies of the type (1.1), as long as f has quadratic growth, and of more general functionals where the surface term has the forḿ Ju ϕ(u + , u − , x, ν)dH n−1 (x) for a suitable surface energy density ϕ : R n × R n × Ω × S n−1 → R. For a discussion of the related problem of density for partition problems we refer to [7].
The restriction of the mentioned results of [9,10,26] to the quadratic energies does not originate from simplicity of presentation, but is instead a consequence of the type of construction used. Indeed, the key idea, first introduced in [9], is to replace the function u by a componentwise linear approximation on a suitably chosen (very fine) cubic grid, and then to remove the cubes which intersect, in a suitable sense, the jump set. The fact that the energy is quadratic permits an explicit integration of the energy density in each cube, and leads to the identification of the continuum energy of the approximation with a discrete energy, which consists of sums of squares of difference quotients along the edges of the grid. In turn, for a suitable choice of the grid this discrete energy approximates the original energy. For nonquadratic expressions the first step, in which one integrates explicitly over a unit cell, breaks down. Estimates are of course still possible, but the result will only hold up to a p-dependent factor, even in the easy case where the functions are smooth to start with. Therefore we use a different strategy, and resort to a piecewise affine interpolation on a suitable refinement of the grid, see discussion in Section 3 below.
Our result permits to replace functions GSBD p (Ω) ∩ L p (Ω; R n ) with much more regular functions in a number of problems related to fracture (see for example [11,[13][14][15]). After the completion of this work, Chambolle and Crismale in [11] have extended our main result to all functions in GSBD p (Ω) by adopting a different technique. We stress that the extra integrability hypothesis that we impose on the relevant function is often not meaningful for problems in fracture mechanics.
Together with the elliptic regularity results for solutions to linear elasticity type systems established in [17], our result has been instrumental for the proof of existence in dimension n = 2 of minimizers for the strong counterpart of the Griffith functional in (1.1), that was presented in [14,15]. More precisely, in [15] it is proved that any local minimizer u of (1.1) has relatively closed jump set, i.e. H 1 (J u ∩ Ω \ J u ) = 0, and it is smooth outside it, namely u ∈ C 1,α (Ω \ J u ; R 2 ) for some α ∈ (0, 1). The equivalence between the weak formulation of the problem as stated in (1.1) and the classical strong form then follows (cf. [14,15] for more details). Such a mild regularity result extends the analogous statement for SBV p functions proved in the celebrated paper [21] by De Giorgi, Carriero and Leaci, corresponding in applications to the (generalized) antiplane shear setting. As already mentioned before, in [14,15] the strong approximation property established in this paper is used to infer such kind of regularity; viceversa the quoted approximation result by Cortesani and Toader [19] uses De Giorgi, Carriero and Leaci's regularity result (in particular by means of [6], Lem. 5.2 by Braides and Chiadó Piat) as a key tool to prove the strong approximation property.
In closing this introduction we mention a complementary approach to the regularity of SBD p and GSBD p functions, which has received a lot of attention in the last years, namely, the proof of rigidity estimates for functions with small jump set. A Korn-Poincaré bound in term of the elastic energy alone was proven for SBD p functions in [12]. An improved estimate, which controls also the gradients, was obtained in the two-dimensional case in [16,23,24].
Finally, we summarize the structure of the paper. Section 2 is devoted to fix the notation for the piecewise affine finite elements and the functional spaces which are involved in our main approximation result, Theorem 3.1, that we shall prove in Section 3.

Notation
One key ingredient of our piecewise affine approximation is the Freudenthal partition of the n-cube [0, 1] n . We say that the vertex (i 1 , . . . , i n ), i k ∈ {0, 1}, precedes the vertex (j 1 , . . . , j n ), j k ∈ {0, 1}, if i k ≤ j k for all k. The convex hull of a chain of n + 1 distinct vertices is a n-simplex. Then the Freudenthal partition S of the n-cube is defined as the set of all n-simplexes obtained through maximal chains of ordered vertices connecting the origin to the vertex (1, . . . , 1) (see Fig. 1). Such partition counts n! simplexes.
Alternatively, for any permutation σ of {1, . . . , n} one defines a simplex S σ as the convex envelope of the points v i := j≤i e σ(j) , i = 0, 1, . . . , n. Explicitly one obtains It is then apparent that S σ consists of the points x ∈ [0, 1] n such that j → x σ(j) is nonincreasing. Therefore the sets S σ have disjoint interiors and cover [0, 1] n . They differ only by a permutation of the components, hence they are congruent and each has volume 1/n!. We use standard notations for the space BV and its subspaces SBV p always referring to the book [2] for details.
An L n -measurable function u : Ω → R n belongs to GSBD(Ω) if there exists a bounded positive Radon measure λ u ∈ M + b (Ω) such that the following condition holds for every ξ ∈ S n−1 : for H n−1 -a.e. y ∈ Ω ξ the function u ξ y (t) = u(y + tξ) · ξ belongs to SBV loc (Ω ξ y ), where Ω ξ y := {t ∈ R : y + tξ ∈ Ω}, and for every Borel set , the aforementioned quantities e(u) and J u are still well-defined, and are respectively integrable and rectifiable in the previous sense. Moreover for every ξ ∈ S n−1 and for H n−1 -a.e. y ∈ Ω ξ we have where (u ξ y ) denotes the absolutely continuous part of the distributional derivative. In analogy to SBD p (Ω), the subspace GSBD p (Ω) includes all functions in GSBD(Ω) satisfying e(u) ∈ L p (Ω; R n×n ) and H n−1 (J u ) < ∞ (cf. [20]).

The main result
Theorem 3.1. Let Ω ⊂ R n be a bounded Lipschitz set and let p > 1. Given u ∈ GSBD p (Ω) ∩ L p (Ω; R n ), there exists a sequence (u j ) ⊂ SBV p ∩ L ∞ (Ω; R n ) such that each J uj is contained in the union S j of a finite number of closed connected pieces of C 1 -hypersurfaces, u j ∈ W 1,∞ (Ω \ S j ; R n ), and the following properties hold: Remark 3.2. The sequence (u j ) in Theorem 3.1 can be constructed in a way that it satisfies in addition These further properties can be obtained by following the arguments in [26] step-by-step, with obvious modifications due to the fact that the proof there is written for p = 2. In this respect, since only a technical effort is required, we focus here on the main difficulties and we prove Theorem 3.1 in the stated form.
As mentioned in Section 1 the proof of Theorem 3.1 follows the general strategy of Chambolle and Iurlano [9,10,26], but uses a different interpolation scheme and a different finite-element grid for the actual construction. Indeed, we first construct a sequence of SBV p ∩ L ∞ -functions converging to a given u ∈ GSBD p (Ω), in a way that the bulk estimate is sharp and the surface estimate is obtained up to a multiplicative factor. Each approximating function is a piecewise linear interpolation outside from a finite number of cubes, where it is set equal to 0. Considering piecewise linear interpolations is essential in order to treat the case of maps in GSBD p (Ω) with p = 2. It is the main difference with the mentioned references [9,10,26] which deal with the quadratic case p = 2. Indeed, in [9,10,26] piecewise polynomial interpolations (of degree equal to the dimension of the space) are employed. Such approximations in dimension higher than 3 if p = 2 would give rise to a multiplicative factor in the bulk estimate (cf. with [9], Lem. A.1), so that the strong approximation property would fail. The piecewise polynomials correspond to a componentwise affine interpolation, that can be done directly on a cubic grid. In the p = 2 case we need to use a piecewise affine interpolation, and therefore need to decompose the domain in simplexes. However, the strategy of [9,10,26] was based on controlling the longitudinal difference quotients along grid segments (the segments joining two vertices of the grid, which are edges of the elements or diagonals of their faces). A natural approach would be to choose an expression for the energy density which uses only these components. In dimension n = 2 this still works, since one can decompose the square [0, 1] 2 into two triangles whose sides have the same orientations (the three orientations being (1, 0), (0, 1) and (1, 1) for both of the triangles). In dimension 3 and higher this is, however, not any more possible, and the energy density will typically not match the geometry of the simplex. Therefore we need to decompose each term of the energy into the components which are "longitudinal" with respect to the edges of the simplexes. We shall denote by A ∈ A ⊂ R n×n sym the "components" of e(u) which enter the energy, and by α A,S j the coefficients of the decomposition of strain direction A in linear combinations of longitudinal strains along the edges simplex S, where j labels the sides of S. The key observation on which the construction in this paper is based is that one can perform this decomposition jointly for the continuous and for the discrete energy.
We now introduce the objects just mentioned in more detail. We fix p ≥ 1 and choose a finite set of matrices A ⊂ R n×n sym , which span R n×n sym and are fixed for the rest of the proof. Let W : R n×n sym → R be defined by where A : B := Tr A T B = ij A ij B ij denotes the Euclidean scalar product on R n×n sym . We denote by D S the set of the edges directions for a simplex S in the Freudenthal partion S . Notice that D S contains n(n + 1) /2 linearly independent vectors and that for any given S the set {e ⊗ e : e ∈ D S } constitutes a basis for R n×n sym . To see this, it suffices to show that if ξ ∈ R n×n sym obeys e · ξe = 0 for all e ∈ D S then ξ = 0. Indeed, the simplex S can be written as the convex envelope of {0, f 1 , . . . , f n }, where (f i ) i=1,...,n is a basis of R n . The set D S is then given by the set of the ±f i 's and the set of all the differences f i − f j 's with i = j. Therefore, if ξ ∈ R n×n sym is such that e · ξe = 0 for all e ∈ D S , we deduce first that f i · ξf i = 0 for all i by taking e = f i , and then f i · ξf j = 0 for all i = j by taking e = f i − f j . Hence, ξ = 0. We stress that D S is a set of differences of vertices of S, not a set of unit vectors.
As in the references mentioned above, the key point is to prove an approximation result that enlarges the jump set by at most a fixed factor. The sharp constant can then be recovered by applying this to the complement of a suitable "large" compact subset of J u .

Theorem 3.4.
Let Ω ⊂ R n be an open bounded set with Lipschitz boundary and let p ≥ 1. Given u ∈ GSBD p (Ω) ∩ L p (Ω; R n ), there exists a sequence (u j ) ⊂ SBV p ∩ L ∞ (Ω; R n ) such that each J uj is contained in the union Σ j of a finite number of (n − 1)-dimensional faces of closed simplexes, u j ∈ W 1,∞ (Ω \ Σ j , R n ), and the following properties hold: W (e(u j )) dx + H n−1 (Σ j ) ≤ˆΩ W (e(u)) dx + c 1 H n−1 (J u ).
In order to prove Theorem 3.4 we need a preliminary lemma, whose proof is entirely similar to Lemma 3.2 of [9] and Lemma 3 of [26] and therefore not repeated here. The given function u is replaced by another GSBD p -function close in energy to u and defined in a larger set.
Fixed y ∈ [0, 1) n and h > 0 small, we consider the translated lattice hy + ξ, with ξ ∈ hZ n . We introduce the tubular neighborhood in the direction −τ of Jû, x − τ ] = {y ∈ R n : [y, y + τ ] ∩ Jû = ∅}, for τ ∈ R n , and the longitudinal difference quotient along the edgeẽ j of S ∈ S S j,h (z) := (û(z + ha j + hẽ j ) −û(z + ha j )) ·ẽ j h|ẽ j | 2 for [z + ha j , z + ha j + hẽ j ] ⊂Ω, and zero elsewhere, where a j and a j +ẽ j are the only two vertices of S whose difference isẽ j . Let us introduce the discrete bulk and surface energies where 1 B denotes the characteristic function of the set B, α A,S j are the coordinates of A in the basis {ν j ⊗ν j : e j ∈ D S } of R n×n sym , whereν j =ẽ j /|ẽ j | andc 1 := 2 n n √ n, the latter choice will be motivated later. Let w y,h be the piecewise affine function obtained interpolatingû on each simplex of the partition. Let us prove that there exist y ∈ [0, 1) n and a subsequence of h ↓ 0 not relabeled, such that (1 ) w y,h −û L p (Ω,R n ) → 0; (2 ) lim h→∞ E y,h 1 (Ω) + E y,h 2 (Ω) ≤ˆΩ W (e(û))dx + c 1 H n−1 (Jû), where c 1 is a constant depending onc 1 and W is the integrand defined in (3.1).
In order to prove (1 ) we observe that for every simplex ξ + hy + hS of the partition with vertices a i , i = 0, . . . , n, there exist n + 1 affine functions f i such that (a i )f i on ξ + hy + hS.
Then integrating on [0, 1) n , we deduce by convexity and Fubini's theorem for h sufficiently small where c takes into account the convexity of the power R t → |t| p and the number of simplexes sharing a certain a i as a vertex. Changing variable z = (x − hy − ξ)/h in the second integral we obtain By dominated convergence theorem the last term vanishes as h ↓ 0, therefore there is a subsequence of h, not relabeled, and a measurable set of full measure E ⊂ [0, 1) n such that for every y ∈ E the convergence in (1 ) holds. Let us prove now (2 ). We integrate again on [0, 1) n and estimate first the bulk part. Changing variable x = hy + ξ and then slicing through Fubini's theorem we obtain where π j is the orthogonal projection on Πν j :=ν ⊥ j . We recall thatν j =ẽ j /|ẽ j | and that the slice is defined as usual byû ν z (s) :=û(z + sν) · ν for z ∈ Π ν . Sinceû ∈ GSBD p (Ω) ∩ L p (Ω; R n ) we haveû ν z ∈ SBV p (Ω ν z ), for H n−1 -a.e. z ∈ Π ν . Observe that 1 J he (z + sν) = 0, for e = |e|ν and z · ν = 0, means z + sν ∈ J he , which is the same as [s, s + h|e|] ∩ (Jû) ν z = ∅. For almost every z, by (2.1) this implies [s, s + h|e|] ∩ Jûν z = ∅. Therefore (3.3) yieldŝ As h ↓ 0 we find by the continuity of the translation, the Lebesgue theorem, and the dominated convergence theorem lim sup Arguing in a similar way for E y,h 2 we obtain having set c 1 :=c 1 #(∪ S∈S D S ) = 2 n−1 n 5 /2 (n + 1)n!, ν e := e/|e|, and having used in the last inequality the slicing formulaˆJû |νû · ν e | dH n−1 =ˆΠ νe #(Jû νe z ) dH n−1 (z) .
By inequalities (3.4) and (3.5) and Fatou's lemma we conclude that there exists y ∈ [0, 1) n and a subsequence of h ↓ 0 not relabeled for convenience such that properties (1 ) and (2 ) hold. In what follows we drop the index y and denote w y,h simply by w h .
We define now a sequence v h as 0 in the cubes Q = ξ + hy + [0, h) n such that Jû crosses an edge of ξ + hy + hS for some S ∈ S , while we set v h := w h otherwise. In the first case we say that the cube is bad, in the second case that it is good. We let Σ h be the union of the faces of the bad cubes. We claim that the constantc 1 in (3.2) can be chosen in a way that for every h sufficiently small As for (2 ), we first notice that H n−1 (J v h ) ≤ 2nh n−1 N h , being N h the number of bad cubes. For every bad cube Q :=ξ + hy + [0, h) n there is at least one pair e ∈ S D S , ξ ∈ hZ n , such that [ξ + hy, ξ + hy + he] ⊂ Q and 1 J he (ξ + hy) = 1. At the same time, the edge [ξ + hy, ξ + hy + he] is shared by at most 2 n−1 cubes. Therefore Recalling that |e| ≤ √ n and definingc 1 := 2 n n √ n we obtain from the definition of E h 2 (Ω) that Let us prove now thatˆΩ Since w h is the affine interpolation ofû on each simplex constituting Q, we have for every pair a, b of vertices of Q, with ν = (a − b)/|a − b|. Thereforê recalling that the difference quotient is defined by and that a j , a j + e j represent the only two vertices of S whose difference isẽ j . Summing on the good cubes Q which intersect Ω we finally obtain (3.7). Property (2 ) then follows by (3.6) and (3.7).
Using a by now standard Besicovitch covering argument we can refine the estimate obtained in Theorem 3.4 reducing the coefficient of the surface term to 1. The idea is to cover the most of the jump set of u with a finite number of pairwise disjoint closed balls, in a way that the jump set in each of them is close to a C 1 hypersurface separating the ball into two components. Then Theorem 3.4 is applied in each component and in the complement of the balls, so that the jump of the resulting function is the union of the C 1 hypersurfaces separating the balls and of the (n − 1)-dimensional faces of closed simplexes obtained by applying Theorem 3.4 (see [9], Thm. 2 or [26], Thm. 6 for more details). Theorem 3.6. Let Ω ⊂ R n be an open bounded set with Lipschitz boundary, p ∈ (1, ∞), and let u ∈ GSBD p (Ω) ∩ L p (Ω; R n ). Then there exists a sequence u j ∈ SBV p ∩ L p (Ω; R n ) such that J uj is contained in the union S j of a finite number of closed connected pieces of C 1 -hypersurfaces, u j ∈ W 1,∞ (Ω \ S j ; R n ), and the following properties hold: (1) u j − u L p (Ω;R n ) → 0; (2) lim sup j→∞ ˆΩ W (e(u j )) dx + H n−1 (S j ) ≤ˆΩ W (e(u)) dx + H n−1 (J u ).
Remark 3.7. We emphasize that the growth hypothesis p > 1 is in fact not needed in the proof of Theorem 3.6. It is only used to deduce conditions (2) and (3) in Theorem 3.1 by means of the GSBD p compactness result in Theorem 11.3 of [20] and strict convexity of W as oulined below.