{'fraction_non_alphanumeric': 0.06198129630379337, 'fraction_numerical': 0.01995870612525809, 'mean_word_length': 3.6206509539842875, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A FETI-DP preconditioner of discontinuous Galerkin method for multiscale problems in high constrast media

In this paper we consider second order elliptic partial differential equations with highly varying (heterogeneous) coefficients on a two-dimensional region. The problems are discretized by a composite finite element (FE) and discontinuous Galerkin (DG) Method. The fine grids are in general nonmatching across the subdomain boundaries, and the subdomain partitioning does not need to resolve the jumps in the coefficient. A FETI-DP preconditioner is proposed and analyzed to solve the resulting linear system. Numerical results are presented to support our theory.


Introduction
We consider the following problem: Find u * ∈ H 1 0 (Ω) such that where Ω ⊂ R 2 is a bounded polygonal domain. We assume that α(x) ≥ α 0 > 0 and α(x) ∈ L ∞ (Ω) which may be discontinuous, while f (x) ∈ L 2 (Ω). The representative examples of the problem (1.1) are subsurface flows in heterogeneous media [18,19] where the heterogeneity varies over a wide range of scales. The aim of this paper is to design and analyze a FETI-DP method for solving such problems based on a composite FE/DG discretization.
Instead of using the full DG method over the whole domain, the composite FE/DG method employs conforming FE methods inside the subdomains, while applies a DG discretization only on the subdomain interfaces to deal with the nonmatching meshes across the interfaces; see [2,5,6,7,11]. The local bilinear forms of the discrete problem are composed of three symmetric terms: the one associated with the energy, the one ensuring consistency and symmetry, and the interior penalty term [25,24] to handle the nonconforming FE spaces across the interfaces; see cf. (2.6)-(2.9). Such discretization allows for nonmatching grids which provides greater flexibility in the choice of mesh partitioning and memory storage. This may be useful particularly when the coefficient varies roughly in one subdomain and mildly in the others.
FETI-DP methods, as well as FETI [15,14] and BDDC [4,20], have been well established as a class of nonoverlapping domain decomposition methods for solving large-scale linear systems. These methods have been used widely for standard continuous FE discretization, and verified to be successful both theoretically and numerically; see [27] and references therein. FETI-DP method was firstly introduced in [13] following by a theoretical analysis provided in [21]. In FETI-DP algorithms, we need a relatively small number of continuity constraints across the interface in each iteration step. The continuity of the solution across the subdomain interfaces is enforced by Lagrange multipliers, while the continuity at the subdomain vertices is enforced directly by assigning unique values. The methods were further improved in [12,17,27] to use the continuity constraints on the averages across the edges on subdomain interfaces. The FETI-DP methods have been developed more recently, and possess several advantages over the one-level FETI method; see cf. [27].
The FETI-DP method was firstly considered for composite FE/DG discretization in [7]. We will follow the same framework as described therein. In [7], the discontinuities of the coefficients are assumed to occur only across the subdomain interfaces. The main purpose of this paper is to extend the methodology to the case where the coefficients are allowed to have large jumps not only across but also along the subdomain interfaces and in the interior of the subdomains. We recall that such problems were investigated in the context of FETI methods in [22,23].
In this paper, we will use the same DG bilinear form as in [11], construct our FETI-DP preconditioner as in [7], and define the components of the scaling matrix as proposed in [22]. For the theoretical aspect, we employ the cut off technique and the generalized discrete Sobolev type inequality, cf. [11], as well as the standard estimates of the edge and vertex functions, cf. [27]. It will be proved that the convergence of the FETI-DP method only weakly depends on the jump of coefficients, i.e., linearly depends on the contrast of the coefficients in the boundary layer. Here we define the boundary layer as the union of fine triangles that touch the subdomain boundaries. We also show that the condition number estimate of the proposed method is quadratic dependence on H/h where H is the subdomain diameter and h is the fine mesh size. This quadratic dependence on H/h can be relaxed to a weaker dependence of H/h(1 + log H/h) 2 under stronger assumptions on the coefficients in the interior of the subdomains.
The remaining part of this paper is organized as follows. In Section 2, we introduce the composite FE/DG formulation of problem (1.1). The FETI-DP method is presented in Section 3. The main results of the paper are given in Section 4 about the analysis of the condition number estimate. Numerical results are provided in Section 5 to confirm the theoretical analysis. In the last section we summarize our findings and discuss certain extensions.
Throughout this paper we denote a Sobolev space of order k by the standard notation H k (Ω) with norm given by · H k (Ω) ; see e.g., [1] for exact definition. For k = 0 we use L 2 (Ω) instead of H 0 (Ω) and write the norm as · L 2 (Ω) . In addition, A B stands for C 1 B ≤ A ≤ C 2 B with positive constants C 1 and C 2 depending only on the shape regularity of the meshes.

DG Discretization
In this section we present the DG formulations of problem (1.1) that will be studied here. Let the domain Ω = ∪ N i=1 Ω i and Ω i be disjoint shape regular polygonal subdomains of diameters H i . Denote the subdomain boundaries by ∂Ω i . For each Ω i , we introduce a shape regular triangulation T h (Ω i ) with the mesh size h i . Note that the resulting triangulation of Ω is in general nonmatching across ∂Ω i .
We assume that the substructures {Ω i } N i=1 form a geometrically conforming partition of Ω, i.e., the intersection ∂Ω i ∩ ∂Ω j (i = j) is either empty, or a common vertex or edge of Ω i and Ω j . Let us denote the common edgeĒ ij =Ē ji := ∂Ω i ∩ ∂Ω j . Although E ij and E ji are geometrically the same object, we will treat them separately since we consider different triangulations onĒ ij ⊂ ∂Ω i and onĒ ji ⊂ ∂Ω j , with the mesh size of h i and h j , respectively. In the text below, we use E ijh and E jih to denote the set of nodal points of the triangulation on E ij and E ji with mesh sizes h i and h j , respectively, andĒ ijh andĒ jih when the endpoints are included. Moreover, the two triangulations T h (Ω i ) and T h (Ω j ) can be merged to obtain a finer but the same mesh onĒ ij and E ji .
We also denote E i∂ := ∂Ω i ∩ ∂Ω when there is an intersection between ∂Ω i and the global boundary ∂Ω. Let us denote by E 0 i the set of indices to refer to the edges E ij , i.e., j of Ω j which has a common edge E ji with Ω i , and by E ∂ i the set of indices to refer to the edges E i∂ . The set of indices of all edges of Ω i is denoted by For simplicity, we assume that the coefficient α(x) ≥ α 0 = 1, which can be fulfilled by scaling (1.1) with 1/ min x α(x). Without loss of generality again, we assume that α(x) is constant over each fine triangle. The analysis will depend on the coefficient in a boundary layer near subdomain boundaries. For each subdomain Ω i , we define the boundary layer Ω h i by i.e., the union of fine triangles in T h (Ω i ) that touch the boundary ∂Ω i . Furthermore, we set Let α i (x) be α(x) restricted to Ω i . We define the harmonic averages along the edges E ij as follows: Note that the functions α ij (x) and h ij are piecewise constant over the edge E ij on the mesh that is obtained by merging the partitions T h (Ω i ) and T h (Ω j ) along this common edge E ij . It is easy to check that Let V h (Ω i ) be the standard finite element space of continuous piecewise linear functions in Ω i . Define The discrete problem obtained by the DG method is of the form: Here each local bilinear form a i (·, ·) is given as the sum of three symmetric bilinear forms: Here ∂ ∂ni denotes the outward normal derivative on ∂Ω i , and δ is a positive penalty parameter. When j ∈ E 0 i , we set l ij = 2, and let α ij and h ij be defined in (2.2). When j ∈ E ∂ i , we set l i∂ = 1, u ∂ = 0, v ∂ = 0, and define α i∂ = α i and h i∂ = h i .
We introduce the bilinear form It is easy to check that d h (·, ·) is symmetric and positive definite, which can induce a broken norm in V h (Ω) by . The next lemma characterizes the equivalence between the bilinear forms a h (·, ·) and d h (·, ·). This equivalence implies the existence and uniqueness of the solution to the discrete problem (2.5), and also allows us to use the bilinear form d h (·, ·) instead of a h (·, ·) for preconditioning. Lemma 2.1. There exists δ 0 > 0 such that for δ ≥ δ 0 and for all u ∈ V h (Ω), we have where γ 0 and γ 1 are positive constants independent of h i , H i , α i (x), and u. For the proof we refer to Lemma 2.1 of [11].

FETI-DP Preconditioner for the Schur Complement Systems
In this section, we will give the formulation of our FETI-DP method using the framework introduced in [27,7].
3.1. Schur Complement Systems and Discrete Harmonic Extensions. Firstly, we borrow the notations from [7]. Let iĒ ji }, i.e., the union of Ω i and theĒ ji ⊂ ∂Ω j with j ∈ E 0 i , and let Then we set We introduce W i (Ω i ) as the FE space of functions defined on the nodal values of Ω i . That is, In the following, we use the same notation to denote both FE functions and their vector representations. The local bilinear form a i (·, ·) in (2.6) is defined over W i (Ω i ) × W i (Ω i ), and the associated stiffness matrix is given by where ·, · denotes the Euclidean inner product associated to the vectors with nodal values. We , where u i,I represents values of u i at interior nodal points on I i and u i,Γ at the nodal points on Γ i . Note that for subdomains Ω i which intersect ∂Ω by edges, the nodal values of W i (Ω i ) on ∂Ω i \Γ i are treated as unknowns and belong to I i . Hence, we can rewrite (3.4) and the matrix A i as where the block rows and columns correspond to the nodal points of I i and Γ i , respectively. The Schur Complement S i of A i , with respect to the nodal points of Γ i , takes the form Note that S i satisfies the energy minimizing property subject to the condition that w i = (w i,I , w i,Γ ) ∈ W i (Ω i ) and w i,Γ = u i,Γ on Γ i . The bilinear form a i (·, ·) is symmetric and nonnegative with respect to W i (Ω i ), see Lemma 2.1. The minimizing function of (3.7) is called the discrete harmonic extension in the sense of a i (·, ·), denoted by H i u i,Γ , and satisfies , the standard discrete harmonic extension in the sense of a i (·, ·), which is defined by Note that the extensions, H i and H i , differ from each other in the sense that H i u i,Γ at the interior nodes I i depends only on the nodal values of u i,Γ on Γ i while H i u i,Γ depends on the nodal values of u i,Γ on Γ i . The next lemma shows the equivalence between H i and H i in the energy form induced by d i (·, ·). This equivalence will allow us to take advantages of all the discrete Sobolev results known for H i discrete harmonic extensions. The fundamental idea of the proof comes from [6], and we still include the proof here for completeness.
Proof. Here and below, for simplicity of presentation, we omit the subscript Γ and denote u i,Γ by u i if there is no confusion.
The left-hand inequality of (3.10) follows from the energy minimizing property of the discrete harmonic extension H i in the sense of a i (·, ·), and the fact that Here we remain to prove the right-hand inequality.
It is easy to verify that Since the bilinear form d i (·, ·) is symmetric and nonnegative, by Cauchy-Schwarz inequality again, we have with arbitrary > 0. Since , using the formulations (2.9) and (3.8), we get which together imply that We proceed the same lines of Lemma 2.1 in [11], and finally obtain where we have used (3.12). Combining (3.14) and (3.15), we have The right-hand side of (3.16) follows by choosing a sufficiently small . Lemma 2.1 and Lemma 3.1 together directly give the following corollary.
Let us introduce the product spaces 2) and (3.4) for the definitions of W i (Ω i ) and W i (Γ i ), and also (3.1) for notation. We also define where S i is given in (3.6).
With a similar decomposition as (3.2), we can partition W i (Γ i ) as can be written as satisfies the continuity condition (3.21). The subspace of W (Ω ) of functions which are continuous on Γ is denoted by W (Ω ); c.f., Definition 3.3 in [7]. Note that there is a one-to-one correspondence between vectors in V h (Ω) and W (Ω ).
Next we define the nodal points associated with the corner variables by We now consider the subspace W (Ω ) ⊂ W (Ω ) and W (Γ ) ⊂ W (Γ ) as the space of functions that are continuous on all the satisfies the continuity condition (3.23). The subspace of W (Ω ) of functions which are continuous on V i is denoted by W (Ω ); c.f., Definition 4.1 in [7].
We can represent u ∈ W (Ω ) as u = (u I , u Π , u ∆ ), where the subscript I refers to the interior degrees of freedom at nodal points I; see (3.1), the Π refers to the primal(Π) variables at the corners V i for all 1 ≤ i ≤ N , and the ∆ refers to the dual(∆) variables at the remaining nodal points on Γ i \V i for all 1 ≤ i ≤ N . Similarly, a vector u ∈ W (Γ ) can be uniquely decomposed as u = (u Π , u ∆ ). Therefore, we can represent W (Γ ) = W Π (Γ ) × W ∆ (Γ ), where W Π (Γ ) and W ∆ (Γ ) refer to the Π− and ∆−degrees of freedom of W (Γ ), respectively.
Let A be the stiffness matrix obtained by restricting the block diagonal matrix A from W (Ω ) to W (Ω ), where A = diag{A 1 , · · · , A N }. Note that the matrix A is no longer block diagonal since there are couplings between primal(Π) variables. Using the decomposition u = (u I , u II , u ∆ ) ∈ W (Ω ), we can partition A as Note that the only coupling across subdomains are through the Π variables where the matrix A is subassembled.
Once the variables of I and Π sets are eliminated, the Schur complement matrix associated with the ∆−variables is obtained of the form Note that S is defined on the vector space W ∆ (Γ ). where the minimum is taken over w = (w I , w Π , w ∆ ) ∈ W (Ω ) with w ∆ = u ∆ .
The proof of the above lemma can be found in Lemma 6.22 of [27] and Lemma 4.2 of [21].
Next we introduce some notations to define the jump matrix B ∆ . The vector space W ∆ (Γ ) can be further decomposed as (3.26) where the local space W i,∆ (Γ i ) includes functions associated with variables at the nodal points of Γ i \V i . Hence, a vector u ∆ ∈ W ∆ (Γ ) can be represented as In order to measure the jump of u ∆ ∈ W ∆ (Γ ) across the ∆−nodes, we introduce the space The jumping matrix B ∆ can be written as where the rectangular matrix B We can reformulate the discrete problem (2.5), on the space of W ∆ (Γ ), as a minimization problem with constraints given by the continuity requirement: Find u * ∆ ∈ W ∆ (Γ ) such that J (u * ∆ ) = min J (v ∆ ), (3.29) where the minimum is taken over v ∆ ∈ W ∆ (Γ ) with constraints B ∆ v ∆ = 0. The objective function where S is defined in (3.25) and where f i is the load vector associated with the subdomain Ω i , and f can be represented as where the entries f i,∆ are defined as Ωi f v i,∆ dx when v i,∆ are the canonical basis functions of W i,∆ (Γ i ).
Note that A and S are both symmetric and positive definite; see also Lemma 3.3. By introducing a set of Lagrange multipliers λ ∈ W ∆ (Γ), to enforce the continuity constraints, we obtain the following saddle point formulation of (3.29): Find u * ∆ ∈ W ∆ (Γ ) and λ * ∈ W ∆ (Γ) such that Once λ * is computed, we can back solve and obtain 3.3. FEIT-DP Preconditioner. We will now define a preconditioner for F in (3.33).
Let us introduce the diagonal scaling matrix D (i) ∆ , which maps W i,∆ (Γ i ) into itself, for all 1 ≤ i ≤ N . Each of the diagonal entries of D (i) ∆ corresponds to one ∆−node, and it is given by the weighted counting function [22] δ † j (x) := where α i is defined in (2.1). Note that one edge is shared by two subdomains. The union of all these functions δ † j (x) provides a partition of unity on all ∆−nodes. We also define An important role will be played by the operator where δ † j (x) is defined in (3.35). Hence, P ∆ preserves jumps in the sense that B ∆ P ∆ = B ∆ , (3.40) which implies that P ∆ is a projection with P 2 ∆ = P ∆ . Define where S i,∆ is the local Schur complement S i , see (3.6), restricted to W i,∆ (Γ i ) from W i (Γ i ), i.e., S i,∆ is obtained from S i by deleting rows and columns associated with the variables at nodal points of V i ⊂ Γ i .
The FETI-DP method is the standard preconditioned conjugate gradient algorithm for solving the preconditioned system

Condition Number Estimate for FETI-DP Preconditioner
The main result of our paper is included in the following theorem, which gives an estimate of the condition number for the preconditioned FETI-DP operator M −1 F . Theorem 4.1. For any λ ∈ W ∆ (Γ), there exists a positive constant C independent of h i , H i , α(x) and λ such that Proof. By the general abstract theory for FETI-DP method, see [21] and Theorem 6.35 of [27], the proof of the lower and upper bound in (4.1) follows by checking Lemma 4.2 and Lemma 4.3 as below, respectively.
Note that P ∆ is a projection which maps W ∆ (Γ ) to itself. By choosing we can easily obtain where we have used (3.40). It follows from Lemma 3.3 that where the first minimum is taken over The next lemma gives us a crucial estimate of the norm of P ∆ .
where β is defined in (4.2) or/and (4.4), and C is a positive constant independent of h i , H i , α(x) and w ∆ .
where the minimum is taken over We can represent w as We define linear functions to approximate w i onĒ ij andĒ ji with j ∈ E 0 i as follows: Note thatŵ ∈ W (Γ ); see (3.21). Therefore, representingŵ = (ŵ Π ,ŵ ∆ ), we have B ∆ŵ∆ = 0. Using the definition of P ∆ (3.37), we have Define v ∈ W (Γ ) to be equal to P ∆ (w ∆ −ŵ ∆ ) at the ∆−nodes, and equal to zero at the Π−nodes. Let us represent (3.19) and (3.20). Using (3.38) and (3.39), it is easy to check that (4.8) and We denote by V h (∂Ω i ) the space of continuous and piecewise linear functions on the local boundaries ∂Ω i . It is obvious that (v i ) i ∈ V h (∂Ω i ). By the definitions of S ∆ and S , (3.41), (3.18), and (4.5), we have With (4.7), to prove (4.6), we need to show that

By Corollary 3.2 it remains to prove
with (4.10) First we consider the term I 2 of (4.10). For j ∈ E ∂ i the proof is trivial due to the specific choices of parameters. For j ∈ E 0 i , it follows from (4.8) and (4.9) that sinceŵ ∈ W (Γ ), and δ † j ∈ (0, 1). Here e is a fine edge on the mesh that is obtained by merging T h (Ω i ) and T h (Ω j ) along E ij . We recall that α ij is constant on each e ⊂ E ij and denoted by α e ij . By summing up, we finally get (4.11) Now we estimate the first term I 1 of (4.10). Here we introduce two semi-norms defined on Γ i as follows: for any u i ∈ W i (Γ i ), and as the discrete harmonic extension in the sense of a i (·, ·). Hence the function • H i u i is the minimizing function of (4.12). Note that H i v i at the interior nodes depends only on the nodal values of v i on Γ i , i.e., the interior of subdomains Ω i . This implies that (4.14) where we have used the second inequality of Lemma 4.1 in [22]. We can write (v i ) i as where I h is the usual Lagrange interpolation operator, and for j ∈ E i the finite element cut-off function θ Eij (x) equals to 1 for all x ∈ E ijh and vanishes on all the other nodes; see Definition 4.2 of [22]. By (3.35), we know that Putting (4.15) into (4.14), we obtain (4.16) As stated in [22] that |I h (θ Eij (w j −ŵ j ) j )| H 1/2 (∂Ωi) |I h (θ Eij (w j −ŵ j ) j )| H 1/2 (∂Ωj ) , since the discrete harmonic extensions from E ij to Ω i and Ω j are equivalent in the corresponding H 1 −seminorms.
Here we refer to Lemma 4.19 of [27] with the two dimensional case, and have (4.17) Sinceŵ i is a convex combination of the values of w i at the end points of E ij , we can employ the generalized discrete Sobolev inequality, c.f. Lemma 4.5 of [22], and obtain (4.18) With the same argument as (4.17), we get Using the inverse inequality, and the L 2 stability of the L 2 projection we have (4.20) since (ŵ j ) i and (ŵ j ) j are linear on E ij and E ji , respectively. Let (w j ) j be the average of (w j ) j on E ji . By (4.42) in [7] we obtain where we have used (4.10) in [22], and the H 1/2 stability of the L 2 projection. Substituting (4.21) into (4.20), we have where we have used (2.2). Putting the above inequality into (4.19), we obtain (4.23) where we used the fact that for all e ⊂ E ij Using the L 2 continuity of the nodal interpolation operator I h , and proceeding the same lines of (4.18), we have and (4.25) Substituting (4.26) and (4.11) into (4.10), we get By the summation of the above inequality for all 1 ≤ i ≤ N and noting that the number of edges of each subdomain can be bounded independently of N , we finally obtain (4.6) with β satisfying (4.2).
Next we consider the special case when the coefficient α(x) in the subdomains Ω i satisfies (4.3) for all 1 ≤ i ≤ N . (4.28) It is well-known that; c.f. [27], since α(x) ≥ α k for all x ∈ Ω k . Using (4.44) in [7], we have (4.30) min(α i , α j )|d ij | 2 where we have used (2.3), and the fact that δ is practically chosen such that δ = O(1). Proceeding with the same lines of (4.29), we can obtain since (w i −ŵ i ) i = 0 at the end points of E ij .
Using the inverse inequality, and the L 2 stability of the L 2 projection we have (4.32) where we have used (4.43) in [7]. Hence, (4.33) This immediately gives With the same arguments as in (4.26) and (4.27), we finally obtain (4.6) with β satisfying (4.4).

Numerical Experiments
Let the domain Ω be a unit square (0, 1) 2 . For the experiments, we partition the domain Ω into 4 × 4 square subdomains. The distribution of coefficients in each example is presented by figures. We use the proposed FETI-DP method for the discontinuous Galerkin formulation (Section 3) of the problem, and iterate with the preconditioned conjugate gradient (PCG) method. The iteration in each test stops whenever the l 2 norm of the residual is reduced by a factor of 10 −6 . The penalty parameter δ is chosen to be 5 in all the experiments.
Example 5.1. In our first example, c.f. left picture of Fig. 1, the coefficient denotes a 'binary' medium with α(x) = α on a square shaped inclusion (shaded region) lying inside one subdomain Ω i at a distance of h from both the horizontal and the vertical edges of ∂Ω i , and α(x) = 1 in the rest of the domain. We study the behavior of the preconditioner as h and α varies, respectively.
It follows from Tab. 1 that the condition numbers are independent of the values of α since the coefficient contrast in the boundary layer is exactly equal to 1. This is consistent with our theoretical results.
Adopting different fine mesh sizes h, we obtain the log-log plot of the condition numbers in terms of H/h for α = 10 6 . The left plot of Fig. 2 shows a dependence worse than linear growth, which is expected to become harder as h goes finer. This confirms the estimate of (4.2) that contains a logarithmic factor besides the linear dependence.
Example 5.2. The distribution of coefficient is shown in the right picture of Fig. 1, with inclusions in two neighbouring subdomains with coefficient values both larger and smaller than in the boundary layers. Similar as the above example, we investigate the dependence of the condition numbers on the mesh ratio H/h. The right plot of Fig. 2 tells us the robustness of the quadratic dependence in the estimate of (4.4). Example 5.3. We employ this example to investigate the dependence of our method on the coefficient variation in the boundary layers. The distribution of the coefficient is depicted in Fig. 3. The coefficient α(x) = α in the edge islands (shaded region), and α(x) = 1 else where.
The numerical results reported in Tab. 3 confirm our theoretical results in Theorem 4.1, i.e., a linear dependence of the condition number on the coefficient variation in the boundary layers. It is worth further investigation to provide techniques to remove this dependence. In [22], the authors used a pointwise weight to define the scaling matrix and finally made the performance of the method completely independent of the coefficient contrast for some special cases. However, there was no theoretical support to explain this robustness and this technique is not valid for the present example either.