Reduced Lagrange multiplier approach for non-matching coupling of mixed-dimensional domains

The definition, analysis and approximation of boundary value problems governed by partial differential equations with incomplete constrains at the boundary is a relevant topic in computational fluid dynamicsHRT96 . The so called do nothing conditions were introduced to accommodate incomplete boundary data. These techniques were applied to geometric multiscale problems, where domains with different dimensionality are coupled together by means of suitable interface conditionsFormaggia2001561 . Since problems of higher dimensionality (PDEs in two or three dimensions) are supplied with interface data of low dimensionality (one dimensional or lumped parameter models), the incomplete information must be supplemented by some modeling assumptions. These issues become more challenging if we consider the coupling of PDEs defined in domains of heterogeneous dimensions that are embedded into each other. These models represent many real physical problems, such as the flow through fractured porous media, the study of fiber-reinforced materials or the modeling of the small circulation (microcirculation) in biological tissues (a representative geometrical configuration of such problems is shown in Figure 1, left panel). In these examples, physical models defined in inclusions characterized by a small dimension can be approximated using models of lower dimensionality, giving rise to coupled problems across multiple spatial scales. As a result, the information that is transferred at the interface can not match, because of the dimensionality gap. We describe these cases as non-matching mixed-dimensional coupled problems.

The main objective of this work is to develop a general mathematical framework for the analysis and the approximation of PDEs coupled by non-matching constraints across different dimensions. Considering the non standard formulation of the coupling conditions, we focus on their enforcement using Lagrange multipliers (LM). The enforcement and approximation of boundary/interface conditions using LM is a central topic in the development of the finite element methodBabuska1973179 ; Bramble ; Pitkaranta ; Boffi2021 among many others. More recently, the LM method has been applied to couple PDEs across interfacesBURMAN20102680 ; Melenk just to mention a few examples of a broad field in the literature.

The novelty of this work with respect to such literature consists in the use of the LM method to couple equations defined on domains with heterogeneous dimensions. For this reason, an essential aspect of the work is to shed light on the interaction between the LM formulation and the restriction/extension operators that govern the transition of PDEs across spatial dimensions. We work in the abstract framework of saddle point problems in Hilbert spaces and their approximation through Mixed Finite ElementsBoffi2013 . An abstract framework based on exterior calculus has recently appeared for the formulation and the approximation of mixed-dimensional (coupled) problemsBoon2021757 ; Boon2022 . We will explore the intersection of this work with ours in future works.

We consider here three main aspects of the problem. First, we analyse under what conditions the stability of the LM formulation is preserved after the application of the dimensional restriction operator. In other words, we will study how PDEs coupled by LM behave when a mixed-dimensional formulation is adopted. Second, we focus on the gap between physically relevant quantities at the interface. At the level of the continuous problem formulation, we introduce an approximation parameter, $N$, that affects the richness of the LM space that matches heterogeneous dimensions (being $N=1$ the fully nonmatching scenario and $N\rightarrow\infty$ the perfectly matched case) and we study how it influences the satisfaction of the boundary constraints, or equivalently the magnitude of the gap between interface unknowns. Third, we analyze the difference between the full-scale formulation and the mixed-dimensional formulation, putting into evidence how the so called dimensionality reduction error varies with the characteristic spatial dimension of the problem and with the approximation parameter, $N$.

With this work, we aim at shedding light on the common mathematical framework that embraces many recent works involving the applications mentioned above. LM formulations for Dirichlet-Neumann type interface conditions for these problems were recently proposedAlzettaHeltai-2020-a ; HeltaiCaiazzo-2019-a ; HeltaiCaiazzoMueller-2021-a . In these works, a three dimensional bulk problem for mechanical deformations is coupled to a one dimensional model for the mechanical behavior of fibers and vessels, respectively. Also, the complementary problem made of thin and slender mechanical structures immersed into a fluid or solid continuum is particularly relevantSteinbrecher20201377 ; Hagmeyer2022 and could be addressed with the proposed tools. A preliminary mathematical study of these problems was recently developedKuchta2021558 ; Kuchta2019375 ; Kuchta2016B962 ; Laurino2019 ; boulakia:hal-03501521 .

After considering the classical Dirichlet-Neumann interface conditions, we also address Robin type transmission conditions. This variant of the problem formulation is particularly significant for multiscale mass transport and fluid mechanics problemsDAngelo2008 ; Cattaneo20141347 ; Kppl2018 ; Koch2020 ; Possenti20213356 . Here, we show how a prototype of these applications can be formulated and analyzed in the framework of perturbed saddle point problems.

We organize the work as follows. In section 2 we introduce the problem formulation for Dirichlet transmission conditions and we recall some fundamental assumptions and results. The reduced Lagrange multiplier formulation of the problem is presented and analyzed in section 3, where we state the well posedness of the problem. The extension to Robin type transmission conditions is also addressed here. We illustrate the behavior of the dimensionality reduction error in Section 4, and in Section 5 we discuss the properties of the restriction and extension operators when the small inclusion can be described as a mapping of a reference domain by means of an isomorphism. Section 6 discusses the particular but very important case of inclusions isomorphic to a cylinder, where we enforce the constraints at the boundary of the inclusion by means of the projection onto a Fourier space truncated to the $N$-th frequency. The larger $N$, the better is the satisfaction of the constraint at the boundary. We introduce the numerical approximation of the problem and its properties in section 7, and finally in section 8 we discuss some numerical experiments in support of this theory.

Consider the following model problem:

where $V$, represented in Figure 1, is a set of (possibly disconnected) immersed inclusions (such as sheets, vasculature networks, or macroparticles, see for example Figure 2, left panel), denoted as $V:=\cup_{i}V_{i}$, and whose boundary is denoted by $\Gamma:=\cup_{i}\Gamma_{i}\equiv\cup_{i}\partial V_{i}=\partial V$.

To develop our approach, we consider a weak formulation of (1a) that uses Lagrange multipliers to impose the given boundary condition on $\Gamma$. In particular, we start by testing (1a) with smooth test functions $v$ defined on the entire domain $\Omega$, and that are zero on $\partial\Omega$

where $n$ denotes the outward normal vector to $V$, and we extend the problem to the entire domain $\Omega$ by adding to it the weak form of $-\Delta u=\tilde{f}$ in $V$, and impose continuity of the function $u$ at the interface $\Gamma$:

Here $\tilde{f}\in L^{2}(\Omega)$ is an arbitrary extension of $f$ in the entire $\Omega$. With a little abuse of notation, from now on we will not distinguish between $\tilde{f}$ and $f$. We denote respectively $w^{+},\,w^{-}$ the outer and inner values of a function $w$ with respect to $\Gamma=\partial V$ along the outer normal direction $n$, and we denote by $[w]=w^{+}-w^{-}$ the jump of $w$ across $\Gamma$. Then equations (1a)-(1b) are equivalent to

where we need an additional condition to impose the value of $u$ on $\Gamma$. The natural way to proceed, is to impose such boundary condition through a Lagrange multiplier, resulting in the following weak problem:

Given $f\in H^{-1}(\Omega)$ and $g\in H^{1/2}(\Gamma)$, find $u\in H^{1}_{0}(\Omega),\,\lambda\in H^{-\frac{1}{2}}(\Gamma)$ such that

With appropriate conditions on the regularity of $\Gamma$, this problem admits a unique solution, and using integration by parts it is straightforward to show that

We use the following notations for Sobolev spaces in $\Omega$ and $\Gamma:=\partial V$:

with which Problem (3) can be represented in operator form as follows:

given $f\in{\mathcal{V}_{\Omega}}^{\prime},\,g\in{\mathcal{Q}_{\Gamma}}^{\prime}$ find $u\in{\mathcal{V}_{\Omega}},\,\lambda\in{\mathcal{Q}_{\Gamma}}$ such that

where

To demonstrate the well-posedness of (3) we first address the main properties of the operators $A$ and $B$, which will be also central in the development of the reduced Lagrange multiplier approach.

By Theorems 1 and 2.3 it follows immediately that Problem 3 is well posed and admits a unique solution $(u,\lambda)$ (see, e.g., Boffi2013 ). By construction, the restriction of $u$ to $\Omega\setminus V\subset\Omega$ satisfies problem (1).

When the inclusion $V$ is slender, i.e., one or more of its characteristic dimensions are small compared to the measure of $\Omega$, it may be convenient to reformulate the problem on a subdomain of $V$ whose intrinsic dimension is smaller than $d$, i.e., a representative surface, curve, or point $\gamma\subseteq V$, see Figure 2. We call $m$ the intrinsic dimension of $\gamma$, and we assume that $|V|=O(\epsilon^{d-m})$, where $\epsilon$ is the smallest characteristic dimension of $V$, and $\epsilon\ll|\Omega|$.

The main idea of such reformulation is rooted on the assumption that the relative measure of the domain $V$ w.r.t. to the measure of $\Omega$ allows one to accurately represent functions in ${\mathcal{Q}_{\Gamma}}^{\prime}=H^{1/2}(\Gamma\equiv\partial V)$ using a Sobolev space ${\mathcal{W}_{\gamma}}$ defined only on $\gamma$ (for example, $H^{s}(\gamma)^{N}$ for some $s\in[0,1]$ and for some $N\geq 1$).

The abstract setting that enables us to perform such dimensionality reduction, is presented in the following assumption.

Assumption 1 is equivalent to asking that the transposed operator $R^{T}:{\mathcal{W}_{\gamma}}\mapsto{\mathcal{Q}_{\Gamma}}$ is a linear and bounding extension operator with trivial kernel, i.e.,

The operator $R^{T}$ defines a closed subspace ${\mathcal{Q}_{\Gamma}^{R}}=\operatorname{Im}(R^{T})\subset{\mathcal{Q}_{\Gamma}}$. The restriction of problem 3 to ${\mathcal{Q}_{\Gamma}^{R}}$ reads:

Given $g\in{\mathcal{Q}_{\Gamma}}^{\prime}$ , find $(u_{R},\lambda_{R})$ in ${\mathcal{V}_{\Omega}}\times{\mathcal{Q}_{\Gamma}^{R}}$ such that

or, equivalently:

Given $g\in{\mathcal{Q}_{\Gamma}}^{\prime}$ , find $(u_{R},\Lambda)$ in ${\mathcal{V}_{\Omega}}\times{\mathcal{W}_{\gamma}}$ such that

where $\lambda_{R}=R^{T}\Lambda$. It is straightforward to see that $\lambda_{R}=[\nabla u_{R}]\cdot n\,\ \textrm{in}\ {\mathcal{Q}_{\Gamma}^{R}}$.

Standard results for mixed approximations (see, for example, Boffi2013 Theorem 4.2.3, under the additional assumption that the operator $A$ is self-adjoint) allows one to conclude that there exists a unique solution $(u_{R},\lambda_{R})=(u_{R},R^{T}\Lambda)$ to Problem (15) (or, equivalently, to Problem (16)). Moreover, the solution of the problem is bounded by:

Let us now define and analyze the error due to the dimensionality reduction induced by the operator $R$, namely $e_{u}=u-u_{R},\ e_{\lambda}=\lambda-R^{T}\Lambda$. Subtracting equation (16) from (5) we obtain,

For any $q_{R}\in{\mathcal{Q}_{\Gamma}^{R}}$ we have the following orthogonality property:

Since ${\mathcal{Q}_{\Gamma}^{R}}$ is a closed subspace of ${\mathcal{Q}_{\Gamma}}$, the space of orthogonal functions to ${\mathcal{Q}_{\Gamma}^{R}}$, named ${\mathcal{Q}_{\Gamma}}^{\perp}$, is such that ${\mathcal{Q}_{\Gamma}}={\mathcal{Q}_{\Gamma}^{R}}\oplus{\mathcal{Q}_{\Gamma}}^{\perp}$. Using the orthogonality property (26) can be written as follows,

Problem (27) becomes: find $e_{u}\in{\mathcal{V}_{\Omega}}$, $e_{\lambda}\in Q^{\perp}$ such that

Then, observing that the infsup stability of $B$ shown in Theorem 2.3 is satisfied by the pair of spaces ${\mathcal{V}_{\Omega}},\,{\mathcal{Q}_{\Gamma}}^{\perp}$, exploiting (Boffi2013, , Theorem 4.2.3) applied to problem (5) we obtain the following bounds:

where $\alpha$ is the coercivity constant of $A$, and $\beta_{B}$ is the inf-sup constant of $B$ (9). We note that all the constants of (29) are independent of $R$, but they may depend on the size of $\Gamma$ through $\beta_{B}$. The dimensionality reduction error and the influence of $R$ only appears through $\|g-Bu_{R}\|_{{\mathcal{Q}_{\Gamma}}^{\prime}}$. Estimates (29) can be regarded as a bound of the dimensionality reduction error with respect to the residual obtained by enforcing the boundary constraint $Bu=g$ on the closed subspace ${\mathcal{Q}_{\Gamma}^{R}}$. Notice, however, that it is not enough to guarantee that $g\in{\mathcal{Q}_{\Gamma}^{R}}^{\prime}$ in order to eliminate the dimensionality reduction error, since the solution $u_{R}$ of the reduced problem (16) does not necessarily coincide with the full order solution $u$, due to possible errors in the representation of the Lagrange multiplier $\lambda$ in the reduced dimensionality setting.

To obtain an estimate that considers also this error, we start form the following equations, obtained again from a combination of (16) and (5):

that is, for any $w\in{\mathcal{W}_{\gamma}}$,

Thanks to the stability of problem (16), namely inequalities (18), we have

Finally, exploiting the continuity of $B^{T}$ and of $R^{T}$, using the triangle inequality and recalling that $w$ is a generic function in ${\mathcal{W}_{\gamma}}$, we obtain,

We note that, in contrast to (29), the constants of (32) depend on the properties of $R$.

We proceed similarly for the case of Robin transmission conditions. By subtracting problem (25) from (21) and exploiting the representation ${\mathcal{Q}_{\Gamma}}={\mathcal{Q}_{\Gamma}}^{R}\oplus{\mathcal{Q}_{\Gamma}}^{\perp}$, we obtain that

From the latter equation, choosing $q\in{\mathcal{Q}_{\Gamma}}^{R}$ we obtain the following relation,

Furthermore, observing that for $\Lambda\in{\mathcal{W}_{\gamma}}$ we have $\langle R^{T}\Lambda,q^{\perp}\rangle=0$, we obtain that the dimensionality reduction error related to the Robin-type transmission conditions satisfy the following problem: find $e_{u}\in{\mathcal{V}_{\Omega}}$, $e_{\lambda}\in{\mathcal{Q}_{\Gamma}}^{\perp}$ such that

which shares the same structure of (21), (with the exception that the second equation is tested on ${\mathcal{Q}_{\Gamma}}^{\perp}$, but the infsup stability of $B$ holds true also in this subspace) and consequently shares the same stability property that is

To conclude this section, we anticipate that we will use both (29) and (32) to derive a priori estimates of the dimensionality reduction error with respect to $\epsilon$. However, we would like to warn the reader that the two results are not equivalently suited for this purpose. The residual type estimates do not require additional regularity to the functions on the right hand side. Conversely, the approximation type estimates leverage on the additional regularity of $\lambda$, that may not be guaranteed, to exploit optimal convergence with respect to $N$.

Let $\gamma\subseteq V$ be the lower dimensional representative domain for $V$ (possibly the union of disjoint connected components). We assume that both $V$ and $\gamma$ are Lipschitz, and define a geometrical projection operator

that maps uniquely each point on $\partial V$ to one point on the lower dimensional $\gamma$. In (36), we indicate with $\mathcal{P}(\Gamma)$ the power set of $\Gamma$ (i.e., the set of all possible subsets of $\Gamma$), and with $\Pi^{-1}$ the preimage of $\Pi$. For $s\in\gamma$, we will also indicate with $|\Pi^{-1}(s)|$ the intrinsic Hausdorff measure of the set $\Pi^{-1}(s)$. In particular the Hausdorff measure $\mathrm{d}{\mathcal{H}(\Pi^{-1})}$ is such that $|\Pi^{-1}(s)|:=\int_{\Pi^{-1}(s)}\mathrm{d}{\mathcal{H}(\Pi^{-1})}$ for all $s\in\gamma$.

We observe that, in principle, for different points $s\in\gamma$, the set $\Pi^{-1}(s)$ could have different intrinsic dimensionality (i.e., it could be a curve, a surface, or a point). We will focus on the situation where the Hausdorff dimensionality of $\Pi^{-1}(s)$ is the same for all $s\in\gamma$. To simplify the notation, we define

and we assume that $|\Pi^{-1}(s)|>0$ and that $|\Pi^{-1}(s)|$ is bounded for almost all $s\in\gamma$, i.e.,

When no confusion can arise, we will also omit from the notation the dependence of the measure $\mathrm{d}{D}$ and of the set $D$ on $s$.

Under the above assumptions, Fubini’s theorem implies that the integral of any absolutely integrable function $f$ over $\Gamma$ can be decomposed as

The projection $\Pi$ induces naturally, though $\Pi^{-1}$, an average operator for absolutely integrable functions $f$ on $\Gamma$ defined, for each $s\in\gamma$, as the average of $f$ over the preimage $D(s)$:

A natural extension operator from $\gamma$ to the whole surface $\Gamma$ can be defined through the projection operator $\Pi$, i.e.,

for any smooth $w$ on $\gamma$. Clearly, the extension operator $\mathcal{E}_{0}$ is the right inverse of the average $\mathcal{A}_{0}$:

since the function $\mathcal{E}_{0}w$ associates to each set $D(s)=\Pi^{-1}(s)\subseteq\Gamma$ the constant value $w(s)$, whose average on $D(s)$ coincides with $w(s)$.

These operators can be generalized to their weighted counterparts by defining

for a given choice of orthogonal weight functions $\varphi_{i}\in H^{1}(\Gamma)\cap C^{0}(\Gamma)$ such that $\varphi_{0}\equiv 1$ (generating the definitions of $\mathcal{A}_{0}$ and $\mathcal{E}_{0}$ above) and such that

where $\delta_{ij}$ is the Kronecker delta, and $c_{i}$ are positive numbers. We now work out some sufficient conditions that allow one to extend the operators above to the Sobolev spaces $H^{s}(\Gamma)$ and $H^{s}(\gamma)$, respectively, for $s\in[-1,1]$.

To simplify the exposition, we assume that $V$ is a single, simply connected, and non self-intersecting inclusion, and we assume that

The tensor product structure of ${\widehat{\Gamma}}$ deriving from Assumption 2 allows one to define a reference projection operator onto ${\widehat{\gamma}}$ by the orthogonal projection on the last $m$ axes in the reference coordinates $\hat{x}$ using the iso-morphism $\Phi$, i.e.: