An Optimally Convergent Coupling Approach for Interface Problems Approximated with Higher–Order Finite Elements

An Optimally Convergent Coupling Approach for Interface Problems Approximated with Higher–Order Finite Elements


In this paper, we present a new numerical method for determining the numerical solution of interface problems to optimal accuracy with respect to the polynomial order of the Lagrangian finite element space on polytopial meshes. We introduce the notion of a virtual interface, and on this virtual interface we enforce that ”extended” interface conditions are satisfied in the sense of a Dirichlet–Neumann coupling. The virtual interface framework serves to bypass geometric variational crimes incurred by the classical finite element method. Further, this approach does not require that the geometric interfaces are spatially matching. Our analysis indicates that this approach is well–posed and optimally convergent in . Numerical experiments indicate that optimal and convergence is achieved.


Interface problems are classical descriptions of transfer phenomena between two physical systems. This transfer of physical quantities between systems is described by the set of equations posed on the interface between two subdomains. Most commonly, the continuity of the solution and their normal fluxes are enforced. Problems of this kind are well–studied and generally known to have unique solutions. The typical challenge of problems of these kind are primarily numerical, i.e., optimality of approximation and efficiency of computation. Common interface problems in coupled dynamical systems are fluid–structure interaction and substrate diffusion between heterogeneous porous media. They also occur in the field of domain decomposition, where artificial interfaces are introduced so that the underlying monolithic problem may be decoupled to reduce the overall computational cost required for determine the numerical solution of the problem.

Two geometric issues arise we attempt to approximate the solution of an interface problem in the classical finite element framework. First, the geometric approximation of the curved interface is nonconforming with respect to domain specified by the continuous problem. As such, the convergence of the numerical solution is limited to second order. The second issue is that different computational models are often merged together to model the dynamics of the coupled system. Because the domains of the individual subproblems are individually meshed in most cases, this may lead to geometrically nonmatching discretizations of the interface given by the continuous problem. In this case, the notion of an interface problem in the continuous setting no longer generalizes to the discrete setting.

To account for the possible geometrically nonmatching interfaces, numerous projection or transfer operators were proposed to move values from one discrete interface to another, often simultaneous preserving conservation properties between the discrete subproblems. However, these methods are well–known to be limited to second order accuracy in both the and norms because of the second order approximation of the polytopial interface to the continuous interface. In addition, it is well–known that unphysical oscillations may occur when utilizing these approaches. Several alternative approaches were developed to address this issue for linear elements. A common theme among these methods is the notion of correcting the discrete energy functional to balance out the excesses and dearths in the energy caused by the gaps and overlaps caused by the meshing of the continuous domains. However, this approach is applicable in the setting of nonoverlapping domain decomposition, and they fail to apply to general interface problems. In addition, obvious challenges arise in this framework when the PDE in question does not have an energy principle.

In this work, we propose a new method, based on the polynomial extension finite element method (PE–FEM), for interface problems that overcomes the geometric challenges that arise due to meshing. The idea behind this approach is to leverage the PE–FEM method by enforcing that the extensions of the subproblems and their extended fluxes are continuous in the weak sense on the curved interface given by the continuous problem. We refer to this interface as a virtual interface because it appears nowhere in the formulation of our method. Instead, the interface conditions are enforced by utilizing pullback operators onto the discrete interfaces. For this reason, we name our method the Method of Virtual Interfaces (MVI). To our knowledge, this is the first method for interface problems that allows for optimal convergence of the numerical solution in both the and norms, with respect to the polynomial order of the underlying approximation space.

The paper is organized as the following. In §2, we will establish the preliminary ideas necessary in this work. The definition of the governing equations is given in §3 and the well–posedness and the optimality of its numerical solution is proven in §4. Computational examples verifying the theoretical claims of this paper are given in §5. And finally, we close our paper with concluding remarks in §6.


Let and denote a bounded, open domain having a boundary . We now subdivide into two separate nonoverlapping partitions, such that and their boundaries, respectively, are also smooth. To ensure that this boundary regularity is satisfied, we require that the interface, is completely embedded in , i.e., . In addition, we will denote and . Finally, on we associate with it the unit outer normal vector and respectively.

In the following, we will make heavy use of the subscript , where . Unless otherwise specified, we will assume that whenever the subscript apperas, e.g. , it implies that “ for ” for narrative convenience. Now let be the set of all triangles whose union is a polygonal approximation of . The mesh size of will then be defined as . In addition, we will denote . We denote and as the set of all edges whose union is a polygonal approximation of and respectively. We remark that, in general, the polygonal approximation of the interface are geometrically nonmatching and therefore we have that typically . We will also define be the polygonal approximations of . Finally, we define as the unit outer normal vector of .

Left: Example of a continuous domain configuration. Right: Example of a discrete domain approximation.
Left: Example of a continuous domain configuration. Right: Example of a discrete domain approximation.
Left: Example of a continuous domain configuration. Right: Example of a discrete domain approximation.
Left: Example of a continuous domain configuration. Right: Example of a discrete domain approximation.

We now consider mappings between the approximate and continuous interfaces. Let be a point in and be a point in . Because of the piecewise linear nature of the approximation to , there inherently exists a –mapping that maps onto . This mapping shall be denoted

In addition, because of the piecewise smooth nature of the mapping , there exists a piecewise smooth composite mapping that maps , where . We denote this mapping as

Let now denote the Euclidean norm. Then, we shall make the following distance assumptions

for some presumably small . We shall assume that for sufficiently small. In most polygonal meshes, this assumption is almost certainly satisfied because it is a well–known result that the distance between a piecewise linear interpolant of a smooth surface is an accurate representation of the smooth surface.

In this paper, we will perform our analysis in the Sobolev spaces defined as the space of functions whose –th weak derivatives are square integrable, where the set can be any subset of or , and is a rational number. Recall for , the norm is defined as

where . In addition, the seminorm is defined as

The norms for the fractional order Sobolev spaces is then defined as the following

where is any subset of . Also, recall that the duality pairing between (resp. ) and (resp. ) is defined by as

The broken Sobolev norms defined over the polygons , and its boundary is defined as

Finally, we define the subspaces

for any .

We now define the discrete spaces of piecewise continuous Lagrange polynomials used in this work. Let be the space of polynomials of order over a simplex . First, let

be the approximation space defined over , and

as the discrete trace space defined over the discrete interfaces. Next, we define the subspaces

We also define the discontinuous finite element space

and the discrete differential operator as follows: , for , for any . otherwise.

We will denote the projection operator onto an arbitrary discrete space as , defined by

On the discrete spaces , and their subspaces, the following inverse inequalities are satisfied

And finally, we will make the approximation theoretic assumption that

is satisfied given that .

Lifting operators will be often used in this work. We will denote as any continuous lifting operator and as its norm. Similarly, we will denote be a discrete bounded operator with as its norm.

We conclude this subsection by establishing that the norms over are equivalent with the following

From the piecewise continuous nature of the mapping, there exists a partition such that on each is smooth. For positive order derivatives, the lemma follows simply by standard coordinate transformation arguments after seeing that the integral can be defined in the piecewise sense over each partition (see [1]). For negative norms, the same coordinate transformation argument as above also hold, except this time we apply this argument to the duality pairing in the definition of the negative norm.

2.1Approximation with polynomial extensions

In this subsection, we construct the underlying approximation ideas that allow us to meaningfully extend functions in from the polygonal subdomain approximation onto the continuous interface . We will use this extension operator to define a set of discrete interface conditions on the continuous interface, so that the numerical solution is optimally accurate with respect to the polynomial order of the finite element space.

Averaged Taylor polynomials For every , let be the element of containing and a family of disjoint star-shaped domains w.r.t. the balls , such that , for , and , where is the symmetric difference of and . We also ask that and . See Figure ? to see how star-shaped domains can be built for a triangular mesh. Using averaged Taylor polynomials we define, for and 1:


is the indicator function for the set . Note that is meaningful only at and is zero otherwise. For any and its pair and we write

For we have If , then in every , is a polynomial of degree , and therefore reproduces exactly in any adjacent to the boundary and it is equivalent to the classical Taylor polynomial. For we can therefore write, for generic :

We take now and

which is well-defined for any and . For convenience, we also define as

Clearly . For vector functions , we introduce the vector operator . We use this notation in particular for gradients of scalar functions (i.e., ). Finally, we shall denote

as the remainder term for the averaged Taylor series. Likewise, we shall denote the vector remainder term as for vector valued functions.

3Statement of the Problem

In this section we will first describe the continuous interface problem whose solution we wish to approximate. We then describe a method, based on the polynomial extension finite element method (PE–FEM), that allows us to approximate the solution to this problem to optimal accuracy in the norm.

3.1The continuous interface problem

Let be an elliptic operator defined as

where . The interface problem, in the strong sense, is then the following:

subject to the interface conditions

where is the forcing function. From the standard theory of partial differential equations, by the smoothness of the data and coefficient , we have that . Because of this, we immediately have that there exists extensions such that . The same holds for and except with the extensions in the spaces and respectively.

Using these extensions, the interface conditions may be rewritten as

by virtue of Taylor’s theorem for averaged Taylor series.

3.2Numerical method

We now describe the method of virtual interfaces. First, let us define the bilinear forms and as


respectively, where is the PE–FEM Neumann bilinear form. Then recall that the PE–FEM Dirichlet approximation in the setting of the Dirichlet subproblem (see PE–FEM paper) is the following: Seek such that

for any data , where is the stabilization parameter. We remark that the inclusion of this parameter is for analytical reasons; in actual computations, this parameter is not required. Likewise, in the setting of the Neumann subproblem, the PE–FEM Neumann approximation is the following: Seek such that

for any data . It was determined that these subproblems are well–posed and stable given that . The coupling method is derived by adding and together and taking


The resulting set of equations is the following: Seek such that

An iterative method We now present a Dirichlet–Neumann method for determining the solution of . It will be instrumental in determining that our set of equations is well–posed through a contraction argument since, cannot be shown to be coercive. In practice, this method can be used in lieu of solving the MVI equations in the monolithic sense. We now define the method: Given an initial guess , seek for such that

where is a relaxation parameter chosen to encourage convergence and be defined as the following

4Numerical Analysis

In this section, we will first prove the uniqueness of the coupled problem , for any pair of data , by means of showing that there exists a unique convergent sequence of functions whose limit solves . Subsequently, we will determine that the coupled solution is, in fact, optimally convergent in . To this end, we first define as and as the coercivity constant for . The results of the analysis presented in this section are heavily dependent on these constants, where we must have be bounded above by some constant small enough. This is the theoretical limitation of this method. In practice, while this ratio does affect the stability of the method, the limitation of this ratio is not as severe as one may suspect.

4.1Uniqueness of

First, let be the difference between subsequent iterations, where again, . In addition, let . Then, it becomes clear that satisfies the following set of equations:

for . As a corollary to the stability bounds derived in (PE–FEM paper Theorems 4 and 5), we present two important lemmas:

From the stability bound for the PE–FEM Dirichlet problem, we have immediately that

after using Lemma 16 of the PE–FEM paper and applying the coordinate transformation norm equivalence inequalities.

From the stability bound for the PE–FEM Neumann problem (c.f. Theorem 4, PE–FEM paper), we immediately have that

after applying the triangle inequality and Proposition ?. We now proceed to bound and . First, using the triangle inequality, we have that

Subsequently, by the projection theorem and the definition of the operator norm, we have by the projection theorem that

after applying standard approximation bounds and Lemma 20 in the PE–FEM paper. Further, again by Lemma 20 in the PE–FEM paper, we have that

after seeing that and applying the inverse inequality.

Additionally, we have that

after applying the projection theorem, the definition of , the continuity of the discrete lifting operator, and the continuity of in the topology. Applying the bounds and in yields the result of this Lemma.

Using these lemmas, we have that

We see that if is sufficiently small, then there exists an