A Simple Finite Element Method for Elliptic Bulk Problems with Embedded Surfaces

In this paper we develop a simple finite element method for simulation of embedded layers of high permeability in a matrix of lower permeability using a basic model of Darcy flow in embedded cracks. The cracks are allowed to cut through the mesh in arbitrary fashion and we take the flow in the crack into account by superposition. The fact that we use continuous elements leads to suboptimal convergence due to the loss of regularity across the crack. We therefore refine the mesh in the vicinity of the crack in order to recover optimal order convergence in terms of the global mesh parameter. The proper degree of refinement is determined based on an a priori error estimate and can thus be performed before the actual finite element computation is started. Numerical examples showing this effect and confirming the theoretical results are provided. The approach is easy to implement and beneficial for rapid assessment of the effect of crack orientation and may for example be used in an optimization loop.


Introduction
New Contributions. In this contribution we consider a basic elliptic problem with an embedded interface with high permeability, which may be used to model the pressure in a medium with cracks or the temperature in composite materials. Our approach is to use a continuous piecewise linear finite element space and simply insert this space into the weak formulation of the continuous problem which consists of a sum of a form on the bulk domain and a form on the interface. Note that the interface cuts through the mesh in an arbitrary way but we avoid using computations on cut elements and instead compensate the lack of regularity across the interface using a mesh which is adapted close to the interface. This approach leads to a scheme which is very easy to implement.
We derive a priori error estimates which shows that the meshsize for elements close to the interface h Γ ∼ h 2 where h is the global mesh parameter used in the bulk mesh. Such a pre-refinement of the mesh leads to optimal order a priori error estimates in terms of the global mesh parameter. Note that no adaptive algorithm is used instead we just split elements that intersect the interface until they are small enough. We start with a quasi uniform mesh and refine to obtain a conforming locally quasi uniform mesh for instance using an edge bisection algorithm.
In forthcoming work we consider schemes using cut elements which does not require adaptive mesh refinement and also works for higher order elements. The method proposed here is however attractive due to its simplicity and may be an interesting alternative in situations where one does need very accurate solutions for instance in the presence of uncertainties or very complicated networks of interfaces, or for optimization purposes.
Earlier Work. The model we use is essentially the one proposed by Capatina et al. [5]. More sophisticated models have been proposed, e.g., in [1,9,10,14], in particular allowing for jumps in the solution across the interfaces. To allow for such jumps, one can either align the mesh with the interfaces, as in, e.g., [11], or use extended finite element techniques, cf. [2,5,7,8]. Our approach, using a continuous approximation, does not allow for jumps, but we shall return to this question in a companion paper.
The approach of superimposing lower dimensional structures independently of the mesh was recently introduced in the context of structural mechanics in [4,6].
Outline. In Section 2 we formulate the model problem, its weak form, and investigate the regularity properties of the solution, in Section 3 we formulate the finite element method, in section 4 we derive error estimates, and in Section 5 we present numerical examples including a study of the convergence and a more applied example with a network of cracks.

Model Problem
Strong Formulation. Let Ω be a convex polygonal domain in R d , with d = 2 or 3. Let Γ be a smooth embedded interface in the interior of Ω without boundary. Then Γ partitions Ω into two subdomains Ω 1 and Ω 2 , where Ω 2 is the domain enclosed by Γ. Let n i be the exterior unit normal to Ω i . See Figure 1.
Consider the problem: find u : Ω → R such that where a| Ω i = a i are given constants, and f ∈ L 2 (Ω), f Γ ∈ L 2 (Γ) are given functions. We also used the notation ∇ Γ = P ∇ for the tangential gradient where P = I − n ⊗ n is the tangent projection. The jump in the primal variable across the interface is defined for ) and that of the normal flux is defined by [n · a∇v] = n 1 · a 1 ∇v 1 + n 2 · a 2 ∇v 2 , where we recall that n 2 = −n 1 on Γ.
Function Spaces. We adopt the usual notation H s (ω) for the Sobolev space of order s on the set ω and we have the special spaces H 1 0 (ω) = {v ∈ H 1 (ω) : v = 0 on ∂ω} and L 2 (ω) = H 0 (ω). For a normed vector space V we let · V denote the norm on V and we use the simplified notation v L 2 (ω) = v ω . We denote the L 2 -scalar product over ω ⊂ R d or ω ⊂ R d−1 by (·, ·) ω . Weak Formulation. Multiplying (2.1) by v ∈ V = H 1 0 (Ω) ∩ H 1 (Γ) and using Green's formula we obtain the weak form where we used the fact that the boundary contributions on ∂Ω vanish due to the boundary condition and then we used (2.2). We thus arrive at the weak formulation: Introducing the energy norm on V , it follows using the Poincaré inequality v Ω ∇v Ω , which holds since v = 0 on ∂Ω, and the trace inequality v Γ v H 1 (Ω 2 ) , that and hence ||| · ||| is a norm on V . The form A is a scalar product on V and by definition A is coercive and continuous on ||| · |||. Therefore it follows from the Lax-Milgram Lemma that there is a unique solution u ∈ V to (2.10).
Regularity Properties. We have the elliptic regularity estimate Observe that by the boundary conditions and the regularity of Using (2.17) we conclude that Thus using elliptic regularity we find that since the right hand side of (2.18) is in L 2 (Γ). Collecting the bounds we obtain the regularity estimate where we note that we have stronger control of u Γ on the subdomains.

The Finite Element Method
To design a finite element method for the problem we use the classical approach restricting the weak formulation (2.10) to a suitably chosen finite dimensional subspace of V . To this end let • T h be a locally quasi uniform conforming mesh on Ω, consisting of shape regular simplices with element size h T and let h = max T ∈T h h T be the global mesh parameter.
• V h be a finite element space consisting of continuous piecewise linear polynomials on T h .
• T h (Γ) denote the set of elements intersected by the interface: The finite element method takes the form: 4 Error Estimates

Preliminaries
• Let ρ be the signed distance function associated with Γ, negative in Ω 1 and positive in Ω 2 . We then have n = ∇ρ where n = n 1 is the unit normal direction exterior to Ω 1 .
• For ζ > 0 let define a tubular neighborhood around Γ by We also have the formula for the so called closest point mapping p : • The tangential gradient is defined by where we recall that P (x) = I − n(x) ⊗ n(x) is the projection onto the tangent plane T x (Γ) to Γ at x.

Interpolation
We introduce the following concepts.
• Let π h : L 2 (Ω) → V h be the Clément interpolant which satisfies the interpolation error estimate where N h (T ) ⊂ T h is the set of all elements which are node neighbors of T .
Note that with this construction we essentially interpolate u e close to Γ and u outside of Γ.
• We consider meshes that are refined in the vicinity of the interface. More precisely, we assume that there are two mesh parameters h Γ and h such that • We chose δ in the definition (4.5) of I h in such a way that which means that χ δ = 1 on N h (T h (Γ)). We note that (4.7) then implies that we may take δ ∼ h Γ in the definition of χ δ .

Remark 4.1
We note that the total number of degrees of freedom N is related to the global mesh parameter as follows Thus we find that for d = 2 we have N ∼ h −2 , which is equivalent to the unrefined mesh, and for d = 3 we have N ∼ h −4 , which is slightly more expensive compared to the unrefined mesh which scales as h −3 .

Lemma 4.1
There is a constant such that Proof. Using the definition of and thus Term I. Using the product rule and the triangle inequality Here we used that by the properties of the extension there holds which we applied with w = ∇ Γ v. Furthermore, we used the bound ∇v Γt (4.20) where Γ t = ρ −1 (t) = {x ∈ R d : ρ(x) = t}, for |t| < δ 0 , followed by the trace inequality where i = 1 for t ∈ [−δ, 0), i = 2 for t ∈ [0, δ], and finally δ ∼ h Γ .
Term II. Using the interpolation error estimate (4.4) we obtain Here we used the estimate that is obtained using similar arguments as in (4.17).
Term III. Using the trace inequality see [12], the interpolation estimate (4.4), the fact (4.7), and finally the stability of the extension we find that Remark 4.2 Alternatively we may use a different extension operator and prove an interpolation estimate which requires less regularity as follows. We include some details for convenience • There is a continuous extension operator We construct v E by first solving the Dirichlet problem ∆v E = 0 in Ω 2 and v E = v on Γ, for which we have the regularity as follows. Term I and II can be estimated in the same way as above. For Term III we have the estimates where at last we used a trace inequality to pass from Γ to Ω 2 .

Error Estimates
Theorem 4.1 The following error estimates hold Proof. (4.44). The proof follows immediately from Galerkin orthogonality and the interpolation error estimate and thus (4.45). For the L 2 estimate we obtain an error representation formula using the dual problem: find φ ∈ V such that where at last we used the elliptic regularity estimate (2.15).

A Convergence Study for a Simple Interface Problem
We consider a problem with f = 0, f Γ = 1, a 1 = a 2 = a Γ = 1 on the domain Ω = (1, e 5/4 ) × (1, e 5/4 ) with a crack at (x 2 + y 2 ) =: r = e. The exact solution to this problem is given as (4 + e) for 1 < r < e, u 2 = 4 − 4e 5 log (r) − 5 4 + 1 for e < r < e 5/4 , and this solution is applied as Dirichlet boundary conditions on ∂Ω, corresponding to a solution depending only on r with u = u 1 = 0 at r = 1 and u = u 2 = 1 at r = e 5/4 . We compare the convergence on a globally refined mesh with a mesh which is locally refined so that h Γ ≤ h 2 at Γ. The convergence is then checked in L 2 norm and H 1 (semi-) norm.
In Figure 2 we show the discrete solution on a given locally refined mesh. We note that optimal convergence is obtained at the cost of locally refining the mesh, Figure 3, whereas a globally refined mesh gives suboptimal convergence in accordance with (4.44) and (4.45), Figure 4.

A More Complex Example with a Bifurcating Crack
In this example we illustrate the modeling capabilities of our approach with application to a more complex problem involving a bifurcating crack.
Model Problem. Let us for simplicity consider a two dimensional problem with a one dimensional crack Γ which can be described as a graph with nodes N = {x i } i∈I N and edges G = {Γ j } j∈I G , where I N , I G are finite index sets, and each Γ j is a curve between two nodes with indexes I N (j). For each i ∈ I N we let I G (i) be the set of indexes corresponding to curves for which x i is an end point. See Figures 5 and 6. The governing equations are given by (2.1)-(2.4) together with two conditions at each of the nodes x i ∈ N , the continuity condition and the Kirchhoff condition We proceed as in the derivation of the weak form in the standard case (2.5)-(2.9). However, when we use Green's formula on Γ we proceed segment by segment as follows where we changed the order of summation and used the Kirchhoff condition (5.2) together with the fact v is continuous to conclude that Thus we conclude that: • The weak formulation is precisely the same in the bifurcating crack case as in the standard case (2.10).
• Since V h ⊂ V the method also takes the same form as in the standard case (3.1) in this more complex situation.
The similar derivation can be performed for a two dimensional bifurcating crack embedded into R 3 , see [13] for further details.
Numerical Example. The crack pattern is modeled using a polygonal chain interpolating higher order curves with each part of the chain of length h/10. The intersection points with element sides are computed and a new polygonal chain containing the old one cut by the intersection points is constructed. In Figure 7 we show the effect on a coarse mesh and on a locally refined mesh. We now compute two different solutions using global refinement and local refinement. We use local refinement at Γ until the smallest meshsize equals that of the globally refined model. In Figure 8 we give the computed solutions using these two approaches. Here a 1 = a 2 = 1 and a Γ = 100, f = f Γ = 0, and we impose, on the domain Ω = (0, 13) × (0, 9.5), u = 1 at x = 0 and u = 0 at x = 13 and homogeneous Neumann boundary conditions at y = 0 and y = 9.5. The corresponding solution with α Γ = 0 is thus a plane.

Concluding Remarks
We suggest a continuous finite element method with superimposed lower-dimensional features modeling interfaces. The effect of these are computed using the higher dimensional basis functions and added to the stiffness matrix so as to yield further "stiffness" to the problem. Due to the fact that we cannot resolve kinks in the normal derivative across the interface we do not obtain optimal convergence orders. We propose a simple adaptive scheme based on an a priori error estimate which guides the choice of optimal local mesh size, to improve the local accuracy, regaining the optimal order of convergence. The resulting scheme is very simple and computationally expedient for many applications such as when optimization of the position of interfaces is of interest.