Well-Posedness and Finite Element Approximation of Mixed Dimensional Partial Differential Equations

We consider a mixed dimensional elliptic partial differential equation posed in a bulk domain with a large number of embedded interfaces. In particular, we study well-posedness of the problem and regularity of the solution. We also propose a fitted finite element approximation and prove an a priori error bound. For the solution of the arising linear system we propose and analyze an iterative method based on subspace decomposition. Finally, we present numerical experiments and achieve rapid convergence using the proposed preconditioner, confirming our theoretical findings.


Introduction
Numerical simulation of diffusion processes in heterogeneous materials is computationally challenging since the data variation needs to be resolved.Thin structures like cracks, fractures, and reinforcements are particularly difficult to handle.It is often advantageous to instead model thin structures as lower dimensional interfaces, embedded in a bulk domain.Mathematically this results in a mixed dimensional partial differential equation (PDE) where the solution has bulk and interface components that are coupled weakly.The aim of the paper is to study well-posedness and regularity of mixed dimensional PDEs, construct and analyze a finite element method for solving the problem, and to develop and implement a preconditioner for efficient numerical solution of the resulting linear system.
Numerical solution to PDEs posed on surfaces is a well established field.An early contribution to finite element approximation of the Laplace-Beltrami equation on surfaces is due to Dziuk in [8], where a finite element approximation is constructed on a polyhedral approximation of the surface.Other approaches include trace based methods as in [6; 5; 18], where functions in the surface are represented by traces of functions in the ambient space of higher dimension.
The review articles [2; 9] survey several methods for solving PDEs on surfaces.Finite element methods for coupled bulk and interface problems are also well studied, see e.g.[10; 15].
Flow in porous media is probably the most prominent application for mixed dimensional PDEs, see [1; 3; 7; 12; 15] and references therein.These methods are fitted, meaning that the full geometry is resolved by the mesh.The bulk is three dimensional and the interfaces are two dimensional surfaces representing fractures.There is also related work treating three dimensional bulk domains with one dimensional embedded structures.Such models include blood vessels embedded in tissue, see [11].
In this paper we consider an elliptic mixed dimensional model problem with a large number of interfaces.The interface model is inspired by the works [3; 17] where a general framework for multidimensional representations of fractured domains is developed, and the work [15] where Robin-type couplings between the bulk and interface are studied.In the bulk and interfaces, we consider Poisson's equation in d and d − 1 dimensions, respectively.In both cases, we use its primal form.We show that the model is well-posed and pay particular attention to problems with a large number of embedded interfaces and how that affects the coercivity bound of the variational formulation of the problem.We further formulate a fitted finite element method and prove a priori error bounds.Additionally, we propose a domain decomposition method based on subspace correction allowing for efficient numerical solution of problems with complex coupled interface structures.In this part we take inspiration from recently developed subspace decomposition algorithms for spatial network models, see [13; 14].We formulate the Schur complement on the union of interfaces and introduce coarse subspaces, using coarse finite element spaces in the bulk interpolated onto the interfaces.Finally, we present numerical examples in two spatial dimensions together with regularity analysis.
Outline: In Section 2 we present the problem formulation and in Section 3 we formulate a fitted finite element method and derive an a priori error bound.Section 4 is devoted to an iterative method for efficient solution of the resulting linear system.Finally, in Section 5 we present numerical results.

Problem formulation
We consider Poisson's equation in the bulk as well as in the interfaces and couple the solutions in bulk and interfaces by Robin conditions.We are interested in the case when there is a large number of interfaces in the domain.Next follows a description of the geometrical setting to accomodate for this partitioning.Then, we present the model problem on weak form and establish its well-posedness.

Mixed dimensional geometry
We consider an open and connected domain Ω of R d (d = 2 or 3) that is partitioned into d + 1 (not necessarily connected) subdomains of different di-mensionalities: where each Ω c (0 ≤ c ≤ d is the codimension) in turn is composed into a finite number of disjoint and connected subdomain segments defined by a set L c , each of which is a planar hypersurface of codimension c.
In the following, we put some requirements on this decomposition of Ω.We define topological subspaces X c of R d for subdomains in codimension c and topological properties such as dense, closed, open, and operators such as closure (cl), boundary (∂) and interior (int) for subdomains of codimension c should be interpreted in the subspace topology of X c .For a set ω, the notation ω denotes the closure of ω in R d .
For codimension c = 0, we let X 0 = Ω, while for c > 0, we define the subspaces X c = X c−1 \ Ω c−1 .We assume that for all ℓ ∈ L c , the subdomain segment Ω c ℓ is open and that their union Ω c is dense.We add the requirement that either This makes it possible to define the adjacency relation E c between subdomain segments i ∈ L c and j ∈ L c+1 by This requirement admits a dense partitioning of a subdomain segment by subdomain segments of one codimension higher.Thus subdomain segment boundary integrals can be expressed as sums of integrals over other subdomain segments.
We focus entirely on the subdomains of codimension 0, 1 and 2. To simplify notation, we introduce I = L 0 , J = L 1 , and K = L 2 .To reduce the complexity of the model, we have chosen to work in a simplified setting: We assume that all subdomains segments Ω c ℓ are planar and polyhedrons, polygonals, lines or points, whichever applies in their respective dimension, and that they have Lipschitz boundary.Lipschitz boundary rules out slit domains and means that, for example, an interface cannot end inside a bulk without connecting to another interface.These assumptions can be relaxed at the cost of a more detailed treatment of the interfaces as is done in [3; 17].To handle curved interfaces, it is possible to consider the planar setting here as an approximation of the curved case.For sufficiently smooth curved interfaces, lift operators on both bulk and interfaces exist that enable an estimate of the error due to such a variational crime, see [10].
In Figure 1, we see an example of a domain fulfilling our requirements.

Model problem
Denote the L 2 -inner product and norm ω , and let L 2 (ω) and H 1 (ω) be the usual Sobolev spaces of L 2 -integrable functions and weak derivatives.If ω is a subdomain segment of codimension greater than Figure 1: An example of a domain Ω fulfilling our requirements.0, then we interpret H 1 (ω) as H 1 ( ω) after a rigid coordinate transformation of ω to a subset ω ⊂ R d−c where the gradient operator, Poincaré inequalites, Green's theorem, and trace theorems in this lower dimensional space R d−c can be applied.For clarity, the symbol ∇ τ is used to denote the gradient operator on subdomain segments of codimension 1.
For functions in the bulk of the domain, we define Occasionally, we refer to components v 0 i of v 0 ∈ V 0 defined as the restriction of v 0 to Ω 0 i .Note that the domain Ω 0 consists of the disconnected union of Ω 0 i for i ∈ I and no continuity between the subdomain segments is assumed.Further, V 0 0 = {v ∈ V 0 : v| ∂Ω = 0} is the space of functions satisfying homogeneous Dirichlet boundary conditions.
For the functions in codimension 1, we define V 1 as follows.The functions in V 1 are strongly enforced to be continuous over common domains of codimension be any function with each v 1 j being the restriction of v 1 to Ω 1 j .Then we define the constrained space We note that V 1 is a Hilbert space, since it is the kernel of the bounded and linear trace operator on the Hilbert space V 1 b .Finally, we define In the bulk domains i ∈ I, let and a right hand side functional F ∈ V * 0 (in the dual space of V 0 ), the model problem can be stated as to find u ∈ V 0 such that for all v ∈ V 0 , ( The bulk and the interfaces are coupled by a Robin-type boundary condition, continuity of the solution between the interface segments is enforced, and a Kirchoff's law type equation for the flux between interface segments applies.Non-homogeneous Dirichlet boundary conditions defined on ∂Ω can be handled by picking any function g ∈ V satisfying the boundary conditions, setting , and obtaining the solution as u + g. Remark 1 (Bulk-interface coupling).For the bulk-interface coupling, we use boundary conditions III as presented in [15] in a porous media flow setting.The interface parameters A j and B j translate to tangential and normal interface permeabilities K 1 τ and K 1 n , and thickness t of the interface as follows: then the weak formulation in (5) corresponds to the strong problem of finding u ∈ V 0 such that with boundary conditions Here n ω is the outward normal defined on the boundary of the subdomains ω.

Well-posedness
To prove existence and uniqueness of solution for the problem (5) we prove that a is coercive and bounded on V 0 and apply the Lax-Milgram theorem.Coercivity cannot be shown directly by a Poincaré or Friedrichs inequality since the subdomain segments are disconnected and not all of them necessarily intersect the domain boundary.We establish coercivity by an iterative procedure, starting to bound norms over subdomain segments far from the domain boundary in norms over subdomain segments and proceed until we reach subdomain segments for which a Friedrichs inequality can be used.
Lemma 1.The bilinear form a is coercive on V 0 , i.e. there is a , and C depends only on the geometry of the subdomain segments {Ω c ℓ } c,ℓ .Proof.In the following, C denotes a constant that is independent of the arbitrary function The value of the constant is generally not tracked between inequalities.Let i 0 ∈ I be a subdomain segment that intersects the domain boundary so that there is a Friedrichs inequality We note that G = (I, J, E 0 ) defines a bipartite undirected graph, where the subdomain segments of codimensionality 0 and 1 (I and J) constitute the vertices and the edges E 0 are the adjacency relations between subdomain segments.It is bipartite because the edges connect vertices in I only with vertices in J and vice versa.Since Ω is connected, i 0 is reachable from all vertices in G. Then it is possible to define a walk through G beginning at an arbitrary i N ∈ I and ending at i 0 .We express the walk as the sequence of vertices it visits by (i N , j N , i N −1 , j N −1 , . . ., j 1 , i 0 ) and require that all vertices in G are visited at least once by the walk.We allow the same vertex to be visited multiple times to guarantee the existence of such a sequence.
We make use of the following inequality, proven in e.g.[4, eq.(5.3.3)],valid for any Further, we note that by Hölder's inequality, we have which together with (7) gives The following trace inequality will be used, valid for For i ∈ I, we have We also have, Since the sequence (i N , j N , . . ., i 0 ) contains all elements in both I and J, we can rewrite the V -norm of v as follows, This procedure is then iterated, using (11) and (12).For the last term v 0 ) .The walk visits the vertices and edges in the graph at most N times, and since we have iterated the steps N times we get where we have added the terms for the edges from E 0 that were not part of the walk.The constant N is the length of the walk and depends on how the subdomain segments are connected with each other and is thus also a constant depending on the geometry of the problem.In the following, we do not track N and let C absorb it.Finally, to obtain the coercivity bound, we use that α ≤ A i , α ≤ A j and β ≤ B j for all i ∈ I and j ∈ J, and obtain from (14) and the definition of a that Lemma 2. The bilinear form a is bounded on V , i.e. there exists a C < ∞ such that for all v ∈ V , and C depends only on the geometry of the subdomain segments {Ω c ℓ } c,ℓ .Proof.Studying the first two sums of a, we note that i∈I The trace inequality from equation (10) gives us that for (i, j) ∈ E 0 , Using this, we can conclude that the third sum of a satisfies where the constant C also depends on the maximum number of interfaces surrounding the subdomains.Theorem 1.Under the above assumptions, equation (5) has a unique solution u ∈ V 0 .
Proof.Since V 0 is a Hilbert space, Lemma 1 and Lemma 2 give coercivity and boundedness of the bilinear form, and the right hand side F is assumed to be a bounded linear functional on V 0 , the Lax-Milgram theorem guarantees existence of a unique solution u ∈ V 0 .

Fitted finite element method
In this section we introduce a fitted finite element discretization of the mixed dimensional model problem and prove an a priori error bound.

Discretization
We let T 1 h,j be a quasi-uniform and shape-regular partition of Ω 1 j , containing closed simplices of diameter no greater than h, with T 1 h = j∈J T 1 h,j .Further, we let T 0 h,i be a corresponding partition of Ω 0 i , with T 0 h = i∈I T 0 h,i and such that the simplices of T 1 h constitute the sides of some of the simplices in Similarly to what we did in Subsection 2.2, for the functions on the interface partition, T 1 h , we define j ) : v 1 piecewise linear functions on T 1 h,j v 1 = 0 on ∂Ω}, and For the functions on the bulk partition, T 0 h , we define and use these definitions to define We note that V 0 h ⊂ V 0 and V 1 h ⊂ V 1 are also Hilbert spaces and the argumentation in Section 2.3 holds for the discretized equation as well.Thus, Lax-Milgram guarantees the existence and uniqueness of a solution to (17).

Interpolation estimates
In this section, we introduce an interpolation operator , where I 0 h and I 1 h are Scott-Zhang [19] interpolation operators onto V 0 h and V 1 h , respectively.For each degree of freedom k (enumerating the Lagrange basis functions {ϕ 0 k } of V 0 h ), we choose an edge (or face if d = 3, etc.) Kk in the support of ϕ 0 k such that the mesh vertex associated with k is contained in Kk .Let ψ 0 k be the L 2 ( Kk )-dual basis, i.e.
and define The interface interpolation operator I 1 h is defined analogously.In the V -norm, we obtain the interpolation error estimate, for v ∈ V 0 , For proofs and further details on this interpolation, we refer to [4, Section 4.8; 19].

A priori error estimate
In order to derive an a priori estimate, we use coercivity and boundedness (Lemma 1 and 2) together with Galerkin orthogonality.For v ∈ V h , it holds that Thus we obtain By choosing v = I h u, and combining ( 18) and ( 19), we can thus state the following theorem.
Theorem 2. (A priori error estimate.)Let u h ∈ V h and u ∈ V be the solutions of (17) and (5).Then .

Iterative solution based on subspace decomposition
Since the bulk domains are only connected through the interfaces it is natural to consider a Schur complement formulation.We propose an iterative solver of the Schur complement equation that is based on subspace decomposition.The subspaces are introduced using an artificial coarse mesh on top the computational domain, see Figure 5 for an illustration and [14; 13] for more details.Given the subdomains we use an additive Schwarz preconditioner to solve for the Schur complement.Given the solution on the interfaces we solve for the decoupled bulk regions, using a direct solver in parallel.

Schur complement
We let {ϕ 0 k } and {ϕ 1 ℓ } be the standard Lagrange bases of V 0 h and V 1 h respectively.Thus the finite element solution u h ∈ V h consists of the two parts We denote by U 0 and U 1 the vectors (U 0 k ) k and (U 1 ℓ ) ℓ , respectively.
Let A and b be the resulting matrix and load vector from ( 17), using the above mentioned basis.Then the system AU = b can be divided into where A 00 and A 11 correspond to the degrees of freedom in the bulk and interfaces, respectively.The submatrices A 01 and A 10 = A ⊤ 01 describe the connections between the interfaces and bulk regions.We form the Schur complement of the submatrix A 00 and note that U 1 solves: The submatrix A 00 is block diagonal which means that linear systems involving this matrix can be solved cheaply in parallel.We therefore focus on how to solve the Schur complement equation ( 21) efficiently.

Preconditioner based on subspace decomposition
The goal is to define a preconditioner for equation ( 21).Inspired by [13; 14], we introduce an artifical quasi-uniform coarse (H > h) finite element mesh T H of the computational domain with corresponding finite element space W H , see Figure 5.We let {φ j } n j=1 be the set of Lagrangian basis functions spanning W H .The mesh T H is independent of the meshes T 0 h and T 1 h and thereby independent of the location of the interfaces.Next, we let I nodal h be the nodal interpolant onto V 1 h , which is the finite element space for the interfaces on which the Schur complement equations are posed.We define the coarse space Furthermore, we define for j = 1, . . ., n, Thereby we have created a subspace decomposition of V 1 h , i.e. any v ∈ V 1 h can be written as a sum v = n j=0 v j with v j ∈ W j .In the remainder of this section, we use matrix representations for functions and operators on the interfaces.The functions {ϕ 1 ℓ } m ℓ=1 are used as basis in V 1 h and the matrices Q j ∈ R m×mj are prolongation matrices mapping from functions in W j to V 1 h .As basis in W j , for j = 0, the functions {φ j } n j=1 are used (i.e.m 0 = n), while for j ≥ 1, we choose the smallest subset of {ϕ 1 ℓ } m ℓ=1 spanning W j , whose size we denote by m j .We now introduce the matrices and form the full preconditioner by the sum, The matrix T is a preconditioner to (21), Note that in order to compute one application of the preconditioner T we need to solve one global coarse scale problem (for j = 0, on the subspace W 0 ) and n independent local problems (for j ≥ 1, on the subspaces W j ).All these solves are done using a direct solver, which means that the technique we propose is semi-iterative.Under some assumptions of the geometry of the interfaces, it is actually possible to show optimal convergence properties of this preconditioner, in the meaning that the convergence rate is independent of the fine scale mesh size h and the coarse scale mesh size H.

Convergence of the preconditioned conjugate gradient (PCG) method
This preconditioner has been thoroughly analyzed for spatial network problems in [13].In order to apply the convergence analysis presented there we need to formulate our Schur complement equation as an equation posed on a graph.We limit ourselves to the case d = 2, so that the interfaces are one-dimensional.We let N be the set of nodes in the interface finite element mesh T 1 h and let E be the set of edges connecting two nodes, where an edge connecting two nodes x and y is represented by an unordered pair {x, y}.This forms the graph G = (N , E).In [13], two operators are used in the analysis: a weighted graph Laplacian L Under Assumption 3.4 in [13], Theorem 3 holds.Moreover, with H ≥ R 0 , Theorem 4.3 in [13] guarantees convergence of the conjugate gradient method for solving Ã11 U 1 = b1 with preconditioner T , defined by equation ( 22).We let U (ℓ) 1 denote the approximate solution after ℓ iteration and let Then the following error bound holds, where κ is the condition number of T Ã11 .Expressed in the V 1 -norm, we have where u 1 .We note that C only depends on upper and lower bounds of the coefficients in equation ( 5), and have from [13] that κ is independent of the mesh sizes h and H. Once the interface component 1 .This problem is easy to solve since A 00 is block diagonal and each block corresponds to a problem in a single bulk subdomain segment.

Numerical examples
We start by investigating the convergence theoretically and numerically for the proposed method for some example problems.We then examine the number of iterations needed for the preconditioned conjugate gradient method to converge, for two different test cases.All numerical examples are two dimensional, d = 2, and posed on the unit square.

Regularity and convergence
The solutions u 0 i on the bulk domains Ω 0 i to equation ( 5) have H 3 /2 -regularity.If the bulk domain is convex it has H 2 -regularity.We formulate this result as a theorem and leave the proof to Appendix A.
In order to ensure that the method converges according to theory, we solve a few problems with several different mesh sizes h.The problems we choose for this analysis are a domain with eight infinite interfaces and a domain with 27 finite interfaces, see Figure 2. Infinite interfaces here means straight lines that pass through the entire domain, leading to convex subdomains, while finite lines leads to non-convex subdomains.For all of these cases, we let A i = 1, A j = 1, B j = 1 and f = e −10 √ 2 )2 on both the interfaces and the bulk areas.The mesh edges are split in two in each refinement of the mesh.The solution using eight interfaces is presented in Figure 4.
The energy norm a(u ref , u ref ) 1/2 of the reference solution u ref is then calculated for each setup and compared to the finite element approximations using different levels of refinement.The initial mesh is fitted to the interfaces and minimal in the sense that interfaces are not refined.We then refine the initial mesh uniformly, using regular refinement, four times.The reference solution is constructed using five uniform refinements.The convergence result is presented in Figure 3.We detect a linear convergence rate, which agrees with the theory for the case of convex subdomains, and a convergence rate that exceeds the predicted one for the case of non-convex subdomains.

Convergence of the iterative solver
We solve equation ( 5) on a unit square with roughly 200 infinite interfaces, using  preconditioner, leads to very poor convergence (around 5 000 -10 000 iterations in the presented examples).
The number of iterations needed to reach convergence using the preconditioner presented in Section 4 for H = 1/8 (the grid in Figure 5) and H = 1/16 are shown in Table 1.We see rapid convergence of the conjugate gradient method given the high condition number for the initial system matrix.We note that a higher number of iterations for the case B j = 100, corresponding to strong coupling, is needed.The convergence rate seems to be independent of, or only depend weakly on H and h, as suggested in the theory.The dense distribution of interfaces and the rapid convergence suggests that Theorem 3 is applicable.It is however difficult to prove this statement since it depends on the exact connectivity and density properties on the network of interfaces.
We also solve the problem on a domain with finite lines, leading to nonconvex subdomains.We do this for roughly 200 lines of length ≤ 0.2.The parameters are the same as in the case of infinite interfaces, so we examine weakly, moderately and strongly coupled cases where A j is either constantly 1 och uniformly distributed between 0.01 and 1.The condition numbers of the full matrix vary in the range 10 8 to 10 13 .The number of iterations using a preconditioner are shown in Table 2. Again, we notice that the preconditioner does a very good job and that the convergence rate is independent or only weakly dependent on the method parameters H and h.We again need a higher number of iterations for the case B j = 100 corresponding to strong coupling.

A Regularity
In this section, we prove that the regularity of a solution to equation ( 5) posed on a two dimensional polygonal Lipschitz domain is H 3 /2 in the bulk regions in general and H 2 in case of convex subdomains.We also show that we have H 2 -regularity on the interface segments.
The strong formulation of the problem posed on With the help of the product rule and the fact that A i ≥ α > 0, this problem can be rewritten as We allow the subdomain Ω 0 i to be non-convex, but assume Ω to be convex.We let Γ k denote the segments of the boundary of Ω 0 i .We let ω k denote the internal angle at corner P k , between Γ k and Γ k+1 .Further, we let D be the set of indices of boundary segments that have Dirichlet boundary conditions, i.e. the segments from ∂Ω.Similarly, we let R be the set of indices of the boundary segments equipped with Robin boundary conditions.We let S be the set of indices k for which Γ k and Γ k+1 both belong to either D or R, and M be the set of indices k for which Γ k and Γ k+1 belong to one each of them.
We have that ω k < π, k ∈ D ∪ M , because of our assumption that Ω is convex.We also have that ω k < 2π for all other k.In Figure 6 we visualize an example of a subdomain Ω 0 i .We are now ready to prove Theorem 4.
: An example of a subdomain Ω 0 i .Angles between different types of boundary conditions have to be smaller than π, while other angles have to be smaller than 2π.
Proof.From Lax-Milgram we have that u 0 i ∈ H 1 (Ω 0 i ) and u 1 j ∈ H 1 (Ω 1 j ).Since u 0 i ∈ H 1 (Ω 0 i ), along with the assumptions that f i ∈ L 2 (Ω 0 i ) and A i is smooth and positive, we have that f i + ∇A i • ∇u 0 i /A i ∈ L 2 (Ω 0 i ).Further, we have that We have that ω k /(2π) / ∈ N for k ∈ S and ω k /(2π) + 1/2 / ∈ N for k ∈ M .Thus, all the requirements of Theorem 2.1 and Proposition 3.1 from [16] are fulfilled for our problem (23) and hence there exist constants α ℓ k depending on f i , A i , B j and u 1 j , such that if λ k,ℓ / ∈ N, which is the only case present in (24).The variables r and θ are the local polar coordinates such that P k = (0, 0), θ = 0 along Γ k+1 and θ = ω k along Γ k .
Isolating one interface segment Ω 1 j , one gets the partial differential equation: which, with Ω ⊂ R 2 , becomes: where x is a parametrisation of the interface segment.Since f j ∈ L 2 (∂Ω 0 i ∩ Ω 1 j ), and B j and A j are smooth, we can conclude that u 1 j ∈ H 2 (Ω 1 j ).

Figure 2 :
Figure 2: Domains in Section 5.1 on which the convergence analysis is performed.

Figure 3 :
Figure 3: The energy norm of the error between approximations u h and the reference solution u ref .

Figure 4 :
Figure 4: The reference solution u ref to the first numerical example with eight interfaces.

Figure 5 :
Figure 5: Domains in Section 5.2 using coarse grid mesh size H = 1/8 for the subspace correction.

Table 1 :
Number of iterations until convergence using the preconditioner, on a domain with approximately 200 infinite interfaces.

Table 2 :
Number of iterations until convergence using the preconditioner, on a domain with approximately 200 finite interfaces.(a)H = 1/8