{'fraction_non_alphanumeric': 0.051375462551275064, 'fraction_numerical': 0.017870973293938785, 'mean_word_length': 3.9032348804500705, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 13, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Analysis and approximation of mixed-dimensional PDEs on 3D-1D domains coupled with Lagrange multipliers

Coupled partial differential equations defined on domains with different dimensionality are usually called mixed dimensional PDEs. We address mixed dimensional PDEs on three-dimensional (3D) and one-dimensional domains, giving rise to a 3D-1D coupled problem. Such problem poses several challenges from the standpoint of existence of solutions and numerical approximation. For the coupling conditions across dimensions, we consider the combination of essential and natural conditions, basically the combination of Dirichlet and Neumann conditions. To ensure a meaningful formulation of such conditions, we use the Lagrange multiplier method, suitably adapted to the mixed dimensional case. The well posedness of the resulting saddle point problem is analyzed. Then, we address the numerical approximation of the problem in the framework of the finite element method. The discretization of the Lagrange multiplier space is the main challenge. Several options are proposed, analyzed and compared, with the purpose to determine a good balance between the mathematical properties of the discrete problem and flexibility of implementation of the numerical scheme. The results are supported by evidence based on numerical experiments.

1. Introduction. In this study we consider coupled partial differential equations on domains with mixed dimensionality, in particular we address the 3D-1D case. The mathematical structure of such problems can be represented by the following formal equations: in Ω, (1.1a) Problem (1.1) can be described as an example of mixed dimensional PDEs. Here, u, u , λ are unknowns, Ω is a bounded domain in R 3 , whereas Λ ⊂ Ω is a 1D structure parameterized in terms of s and d s is the derivative with respect to s. The term λδ Λ is a Dirac measure such that Ω λ(x)δ Λ v(x) dx = Λ λ(t)v(t) dt for a continuous function v and T Λ : Ω → Λ is a suitable restriction operator from 3D to 1D.
Using models based on mixed dimensional PDEs is motivated by the fact that many problems in geo-and biophysics are characterized by slender cylindrical structures coupled to a larger 3D body, where the characteristic transverse length scale of the slender structure is many orders of magnitude smaller than the longitudinal length. For example, in geophysical applications the radii of wells are often of the order of 10 cm while the length may be several kilometers [25,26]. Similarly, in applications involving the blood flow and oxygen transport of the micro-circulation the capillary radius is a few microns, while simulations are often performed on mm to cm scale, with thousands of vessels [3,13,16,29]. Finally, in neuro-science applications a neuron has width of a few microns, while its length is much longer. For example, an axon of a motor neurons may be as long as a meter. Hence, at least 4 orders of magnitude in difference in transverse and longitudinal direction is common in both geo-physics, bio-mechanics and neuro-science. Meshes dictated by resolving the transverse length scale in 3D would then possibly lead to the order of 10 12 degrees of freedom. Even if adaptive and strongly anisotropic meshes are allowed for, the computations quickly become demanding if many slender structures and their interactions are under study.
From a mathematical standpoint, the challenge involved in problem (1.1) is that neither T Λ nor δ Λ are well defined. That is, without extra regularity, solutions of elliptic PDEs only have well defined traces of co-dimension one. Here, T Λ is of codimension two, mapping functions defined on a domain in 3D to functions defined along a 1D curve. The challenge of coupling PDEs on domains with high dimensionality gap has recently attracted the attention of many researchers. The sequence of works by D'Angelo, [9,10,11] have remedied the well-posedness by weakening the solution concept. The approach naturally leads to non-symmetric formulations. An alternative approach is to decompose the solution into smooth and non-smooth components, where the non-smooth component may be represented in terms of Green's functions, and then consider the well-posedness of the smooth component [15]. The numerical approximation of such equations has been also studied in a series of works. The consistent derivation of numerical approximation schemes for PDEs in mixed dimension is addressed in [5]. Concerning approximability, elliptic equations with Dirac sources represent an effective prototype case that has been addressed in [4,17,19], where the optimal a-priori error estimates for the finite element approximation are derived. Furthermore, the interplay between the mathematical structure of the problem and solvers, as well as preconditioners for its discretization has been studied in details in [21] for the solution of 1D differential equations embedded in 2D, and more recently extended to the 3D-1D case in [20].
Stemming from this literature, in this work we adopt and analyze a different approach, closely related to [18,22]. That is, we exploit the fact that Λ is not strictly a 1D curve, but rather a very thin 3D structure with a cross-sectional area far below from what can be resolved. With this additional assumption, we show that robustness with respect to the cross-sectional area can be restored. The major novelty with respect to the previous works is that we address essential type coupling conditions, namely Dirichlet-Neumann conditions, rather than natural type ones, such as the Neumann-Robin or Robin-Robin cases. These coupling conditions pose additional difficulties as the conditions are not a natural part of the weak formulation of the problem. We overcome this difficulty by resorting to a weak formulation of the Dirichlet-Neumann coupling conditions across dimensions by using Lagrange multipliers.
Although the focus of the present work is mostly about analysis and approximation of the proposed approach, we stress that it aims to build the mathematical foundations to tackle various applications involving 3D-1D mixed dimensional PDEs, such as FSI of slender bodies [24], microcirculation and lymphatics [27,30], subsurface flow models with wells [8] and the electrical activity of neurons.

Preliminaries.
Let the domain Ω ⊂ R 3 be convex and composed of two parts, Ω and Ω ⊕ := Ω \ Ω . Let Ω be a generalized cylinder, c.f. [14], that is; the swept volume of a two dimensional set, ∂D, moved along a curve, Λ, in the three-dimensional domain, Ω, see for For simplicity, let us assume that λ (s) = 1 such that the arc-length and the coordinate s coincide. Further, let D(s) = [x(r, t), y(r, t)] : (0, R(s)) × (0, T (s)) → R 2 be a parametrization of the cross section and Γ be the lateral surface of Ω , i.e. Γ = { ∂D(s) | s ∈ Γ}, while the upper and lower faces of Ω belong to ∂Ω. We assume that Ω crosses Ω from side to side. Finally, |·| denotes the Lebesgue measure of a set, e.g. |D(s)| is the cross-sectional area of the cylinder. In general, |D(s)| must be strictly positive and bounded. According to the geometrical setting, we will denote with v, v ⊕ , v , v , functions defined on Ω, Ω ⊕ , Ω , Λ, respectively.
Let D be a generic regular bounded domain in R 3 and X be a Hilbert space defined on D. Then (·, ·) X and · X denote the inner product and norm of X, respectively. The duality pairing between the X and its dual X * is denoted as ·, · . Let (·, ·) L 2 (D) , (·, ·) D or simply (·, ·) be the L 2 (D) inner product on D. We use the standard notation H q (D) to denote the Sobolev space of functions on D with all derivatives up to the order q in L 2 (D). The corresponding norm is · H q (D) and the seminorm is | · | H q (D) . The space H q 0 (D) represents the closure in H q (D) of smooth functions with compact support in D.
Let Σ be a Lipschitz co-dimension one subset of D. We denote with T Σ : H q (D) → H q− 1 2 (Σ) the trace operator from D to Σ. The space of functions in H 1 2 (Σ) with continuous extension by zero outside Σ is denoted H 1 2 00 (Σ) and we remark that 00 (Σ)) * We will frequently use inner products and norms that are weighted. The L 2 and H 1 inner products weighted by a scalar function w, which is strictly positive and bounded almost everywhere, are defined as follows whereas a weighted fractional space H s (Γ; w) is defined in terms of the interpolation of the corresponding weighted spaces. For the norm of such spaces, we introduce the Riesz map S such that for u, v ∈ H 1 (Γ) we have Then S is a compact self-adjoint operator. Assuming that {λ k } k is the set of eigenvalues, {φ k } k the set of eigenvectors of S orthonormal with respect to the inner product 3 (·, ·) L 2 (Γ),w and u ∈ H 1 (Γ) can be expressed as u = k c k φ k then The space H s 00 (Γ; w) is defined analogously, but with S above defined in terms of the H 1 0 inner product. Owing to the positivity and boundedness of w the weighted spaces equal the corresponding non-weighted spaces as sets, but their norms are different.
Central in our analysis are the transverse averages w, w defined as, where dω, dσ, dγ are the generic volume, surface and curvilinear Lebesgue measures. Clearly, Analogously, for functions defined on Λ and Ω respectively, we let d s and ∂ s be the ordinary and partial derivative with respect to the arclength. The operator obtained from a combination of the average operator (·) with the trace on Γ will be denoted with T Λ = (·) • T Γ , as it maps functions on Ω to functions on Λ. Further, let the extension operator E Γ : H 00 (Γ) be defined such that (E Γ v )(x) = v (s), for any x ∈ ∂D(s). Then, the following identity shows that the transversal uniform extension operator is the inverse of the transversal average, With the above notation we are now able to formulate the precise weak formulations of the problems (1.1), which we will call the Problem 3D-1D-1D. The problem reads: given f ∈ L 2 (Ω), g ∈ L 2 (Ω ), q ∈ H 1 2 In addition to the 3D-1D-1D problem we will also consider an intermediate problem where the 3D and 1D problems are coupled at an intermediate 2D surface encapsulating the 1D structure. The strong form is: The corresponding weak formulation of (2.3), referred to as the Problem 3D-1D-2D, reads: given f ∈ L 2 (Ω), g ∈ L 2 (Ω ), q ∈ H 1 2 We conclude this section with the analysis of a fundamental property for the problem formulation that we will address, namely, the characterization of the regularity of the operator T Λ . More precisely we aim to show that T Λ : . This is a consequence of the following lemma.  .
The proof of 2.1 is based on the representation of fractional norms in terms of the spectrum of the Laplace operator and subsequent standard arguments in harmonic analysis. The full proof is reported in the appendix for the sake of clarity.
Under the geometric assumptions stated above for Ω, Γ, Λ, Lemma 2.1 implies the following result.  . Furthermore, from the above Corollary, it is clear that T Λ : H 1 0 (Ω) → H 1 2 00 (Λ). 3. Saddle-point problem analysis. Let a : X × X → R and b : X × Q → R be bilinear forms. Let us consider a general saddle point problem of the form: find u ∈ X, λ ∈ Q s.t.
The Brezzi conditions [6] ensure that the problem (3.1) is well-posed. For our purpose here, we use the following relaxed version of the Brezzi conditions: Here α, β, C, and D are positive numbers.
Here, the coercivity condition (3.2) applies to X, which is a relaxation of Brezzi's original conditions.
We prove that the conditions of Theorem 3.1 are fulfilled choosing Proof. We need to establish the four Brezzi conditions. The bilinear form a(· , ·) is clearly bounded and coercive since for .
To show the inf-sup condition, we will employ a lifting operator, H Ω , from H 1/2 00 (Γ) to H 1 (Ω). In [28] it is established that extension operators for domains having small geometric details (see also [22] for a direct application to this case) there exists a lifting operator We notice that the trace operator is surjective from H 1 0 (Ω) to H H 1 2 00 (Γ), we can find v = H Ω ξ. Using the stability of the harmonic extension we obtain where in the last inequality we exploited the fact that We prove that the hypotesis of Theorem 3.1 are fulfilled with the following spaces Let us consider X equipped again with the norm |||[·, ·]||| and Q equipped with the norm · . Then, we have the following lemmas.
The inf-sup condition holds. We choose v = 0 and obtain For any q ∈ H 1 2 00 (Λ), we consider its uniform extension to Γ named as E Γ q and then we consider the harmonic extension Moreover, using Lemma 2.1 we obtain Therefore, we conclude the proof with the following inequalities, 4. Finite element approximation. In this section we consider the discretization of the Problems 3D-1D-2D and 3D-1D-1D by means of the finite element method. We address two main challenges; first we aim to identify a suitable approximation space for the Lagrange multiplier and to analyze the stability of the discrete saddle point problem; second we aim to derive a stable discretization method that uses independent computational meshes for Ω and Λ, not necessarily conforming to Γ. Let us introduce a shape-regular triangulation T Ω h of Ω and an admissible partition T Λ h of Λ. We analyze two different cases: the conforming case, where compatibility constraints are satisfied by T Ω h and T Λ h with respect to Γ and consequently h = h; and the non conforming case, where it is possible to choose T Ω h and T Λ h arbitrarily. The discrete equivalent of (3.1) reads as finding where with little abuse of notation we use h as the sub-index for all the discretization spaces. This discrete problem is well-posed if the (3.2)-(3.5) conditions applies to X h and Q h . Since X h ⊂ X and Q h ⊂ Q, (3.2)-(3.4) follow immediately and only the infsup condition needs consideration. We summarize this proposition in the Corollary below.
Corollary 4.1. [12, Theorem 2.42] Let X h ⊂ X, Q h ⊂ Q, a(·, ·) and b(·, ·) satisfy the conditions (3.2)-(3.4) then the problem (4.1) is well-posed if the discrete counterpart of (3.5) is satisfied, i.e. there exists a constant β h > 0 such that 4.1. Analysis of the case where T Ω h conforms to T Λ h and to Γ. As conformity conditions between T Ω h , T Λ h and Γ, we require that the intersection of T Ω h and Γ is made of entire faces of elements K ∈ T Ω h . Furthermore, we also set a restriction between T Ω h and T Λ h . We assume that Λ is a piecewise linear manifold. We want that the intersection of Γ with any orthogonal plane to Λ that crosses Λ at the internal nodes of T Λ h , consists of entire edges of T Ω h . As a result of the latter condition we have h = h. For this reason, we denote as T Λ h the mesh on Λ from now on throughout this section.
, with k > 0, the conforming finite element space of continuous piecewise polynomials of degree k defined on Ω satisfying homogeneous Dirichlet conditions on the boundary and by X k h,0 (Λ) ⊂ H 1 0 (Λ) the space of continuous piecewise polynomials of degree k defined on Λ, satisfying homogeneous Dirichlet conditions on Λ ∩ ∂Ω. The space Q h must be suitably chosen such that (4.2) holds. Let Q h be the trace space of X k h,0 (Ω), namely the space of continuous piecewise polynomials of degree k defined on Γ which satisfy homogeneous Dirichlet conditions on ∂Ω. As a result, In what follows, we analyze the well-posedness of the discrete problem. From now on, C denotes a generic constant independent of the mesh size.
Proof. We show that P h is continuous on L 2 (Γ) and on H 1 0 (Γ) following [12, Section 1.6.3]. Then, Lemma 4.2 can be proved by interpolation between spaces, since H 1 2 00 (Γ) can be seen as the interpolation space between L 2 (Γ) and H 1 0 (Γ). For the L 2continuity, we exploit the fact that, from the definition of Let us now consider v ∈ H 1 0 (Γ). The Scott-Zhang interpolation operator SZ h from H 1 0 (Γ) to Q h satisfies the following inequalities, Therefore, using (4.4), (4.5), the L 2 stability of P h and the inverse inequality, we obtain, from which we obtain the continuity in H 1 0 (Γ). Lemma 4.3. There exists a constant γ > 0 such that for any Proof. From the continuous case, in particular from (3.6), we have and by the trace inequality [1, 7.56]), we obtain sup v∈H 1 0 (Ω) .
Using Lemma 4.2 we obtain, . Theorem 4.4 (Discrete inf-sup). The inequality (4.2) holds true, namely there exists a positive constant β h,1 such that, Proof. As in the continuos case, we choose v h = 0 and we have Using Lemma 4.3 and the boundedness of the harmonic extension operator H Ω from H 1 2 00 (Γ) to H 1 0 (Ω) introduced in the previous section, we have . .
Now we use the conformity of T Ω h to the interface Γ to guarantee that the operator T Γ R h H Ω coincides with the identity on the space Q h . Then, for any q h ∈ Q h we have q h = T Γ R h H Ω q h and owing to this property we obtain the following inequality, which proves the condition, with β h,1 = γ H Ω In this case, we use the same spaces X k h,0 (Ω), X k h,0 (Λ) defined previously. For the multiplier space we choose Q h = X k h,0 (Λ), therefore we impose homogeneous Dirichlet boundary condition on Λ ∩ ∂Ω also for the Lagrange multiplier. We aim to find Below we address the well-posedness of the 3D-1D-1D discrete problem with this alternative choice of multiplier space.
Proof. Again, we choose v h = 0, so that the proof reduces to showing that there exists β h,2 such that For any w ∈ H (Ω) to X k h,0 (Ω). Then, we exploit the conformity of the meshes on Ω, Γ and Λ and the fact that the operator T Λ R h H Ω E Γ coincides with the identity if applied to functions in Q h . As a result, from Lemma 4.6, we obtain the following inequality that concludes the proof with β h,2 = γ H Ω −1 C −1 R .

Analysis of the case where T Ω
h and T Λ h do not conform to Γ. We analyze now the case in which the elements of the 3D mesh T Ω h do not conform with the surface Γ nor with Λ. As the 3D-1D-1D formulation is more suitable for this purpose, we solely focus on the analysis of the discrete version of Problem 3D-1D-1D.
h,0 (Ω) be the approximation of the 3D problem and let u h ∈ X 1 h,0 (Λ) the one of the 1D problem. In contrast to the conforming case, here we limit the analysis to the case of piecewise-linear finite elements. With little abuse of notation, we use the sub-index h for the product space X h = X 1 h,0 (Ω) × X 1 h,0 (Λ). Concerning the multiplier space, let G h = {K ∈ T Ω h : K ∩ Λ = ∅}, be the set of the 3D elements that intersect Λ. Then we define Q h = {λ h : λ h ∈ P 0 (K) ∀K ∈ G h }. We notice that the multiplier functions are defined on the 3D elements. Again with a little abuse of notation, we denote with Q h also the restriction to Λ of the space of piecewise constant functions defined in 3D. As a result, we have Q h ⊂ L 2 (Λ) ⊂ H − 1 2 (Λ). However, with this choice of multipliers the problem is not inf-sup stable, therefore the idea is to add a stabilization term s(λ h , µ h ) : Q h × Q h → R to (4.7a) following the approach introduced in [7]. The objective of this section is to analyze the stabilized version of the 3D-1D-1D problem: find [u h , u h ] ∈ X h and λ h ∈ Q h such that The idea of the stabilization strategy proposed in [7] is to identify a new multiplier space Q H , which is never implemented in practice, such that inf-sup stability with X h holds true. Then, the stabilization operator is designed to control the distance between Q h and Q H through the following inequality being π H a suitable projection operator Q h → Q H . Applying the results obtained in [7], the well posedness of problem (4.9) is governed by the following lemma.
For the proof of this result we refer the reader to Lemma 2.3 of [7]. In the remainder of this section, we show how to find a multiplier space Q H and a stabilization operator s h such that all the assumptions of Lemma 4.8 are satisfied.
The first step consists of showing that there exists a discrete space Q H that satisfies the first assumption of Lemma 4.8. We recall that in the case of Problem The construction of the inf-sup stable space Q H is based on macro elements of diameter H, where H is sufficiently large. In particular, we assume that there exists positive constants c h and c H such that c h h ≤ H ≤ c −1 H h. The space is constructed assembling the 3D elements of G h into macro patches ω j such that H ≤ |ω j ∩ Λ| ≤ cH with H = min j |ω j ∩ Λ| and c ≥ 1. Let M j be the number of elements of the patch ω j , namely, ω j = ∪ Mj i=0 K i , where K i ∈ G h . We assume that M j is uniformly bounded in j by some M ∈ N and that the interiors of the patches ω j are disjoint. We define Q H as the space of piecewise-constant functions on the patches, namely Q H = µ H : µ H ∈ P 0 (ω j ) ∀j . As previously pointed out for Q h , we denote with Q H also the restriction of the multiplier space to Λ, namely say Q H ⊂ L 2 (Λ) ⊂ H − 1 2 (Λ). Moreover, we associate to each patch ω j a shape-regular extended patch (using the classical definition of shape-regularity, see for example [12]), still denoted by ω j for notational simplicity, which is built adding to ω j a sufficient number of elements of T Ω h and we assume that the interiors of the new extended patches ω j are still disjoint (see Figure 4.1). The extended patches ω j are built such that they fulfill the conditions meas( where Γ ωj ∩Λ is the portion of Γ with centerline ω j ∩ Λ. The latter assumption is required to ensure that the intersection of Γ ωj ∩Λ and ω j is not too small and it will be needed later on to prove the inf-sup stability of the space Q H in Lemma 4.9. A representation of this construction in the simple case in which ω j is composed just by one tetrahedron is shown in Figure 4.1. Thanks to the shape regularity of these extended patches, the following discrete trace inequality holds true for any function v ∈ H 1 (ω j ), (Right) the intersection between Γ ω j ∩Λ and ω j . Here for simplicity ω j is represented as a single tetrahedron but actually it is a collection of tetrahedra as shown in left panel.
Moreover, ∀u h ∈ X 1 h,0 (Ω) we have the following average inequality, which is a consequence of the definition of T Λ , Jensen inequality, and the fact that the patches are disjoint We are now ready to prove that the space Q H is inf-sup stable.
Lemma 4.9. The space Q H is inf-sup stable, namely there exists β > 0 such that Proof. We choose v h = 0 and we prove that .
Proving the last inequality is equivalent to finding the Fortin operator π F : H 1 0 (Ω) → X 1 h,0 (Ω), such that 14) We define We notice that supp(T Γ ϕ j ) ⊂ Γ ωj ∩Λ ∩ ω j ensures that T Λ ϕ j ⊂ ω j ∩ Λ. Therefore, since the interiors of ω j ∩ Λ are disjoint and ϕ j = 0 on ∂ω j , the functions T Λ ϕ j ∀j have all disjoint supports. Provided H is sufficiently larger that h, the functions ϕ j and their traces T Γ ϕ j have a sufficiently large support thanks to the fact that meas(ω j ) = O(H 3 ) and diam(Γ ωj ∩Λ ∩ ω j ) = O(H). Owing to these properties it is possible to satisfy (4.16). Then, by construction, Concerning the continuity of π F , we exploit the assumptions that the interiors of ω j are disjoint, supp(ϕ j ) ⊂ ω j and the H 1 -stability of I h to show that For the second term, using that ∇ϕ j L 2 (ωj ) = O(1), ωj ∩Λ |∂D|T Λ ϕ j = O(H) and that |ω j ∩ Λ| ≤ cH, exploiting Jensen's average inequality (4.13) and trace inequality (4.12), and finally applying the approximation properties of I h , the following upper bound holds true (where all the constants have been condensed into C), that is the H 1 -stability of π F . We notice that the constant in the inequality (4.15) is independent of how Λ cuts the elements of the mesh T Ω h . For the second assumption of Lemma 4.8, we recall that b(v h , µ h ) is continuous with respect to the norms |||v h |||, µ h L 2 (Λ) . Using Lemma 4.9, and in particular the existence of a Fortin projector, there exists a constant β such that (the proof is analogous to the one of Lemma 2.1 in [7]) We define π H = j π j H : L 2 (Λ) → Q H , where π j H is the operator Since ∪ j ω j ∩Λ = Λ and ω j ∩Λ are not overlapping, we obtain that π H is an orthogonal projection, namely (w − π H w, π H w) = 0. Moreover, for any w ∈ L 2 (Λ) the following Poincar inequality holds true, see for example [12,Corollary B.65], We consider the following stabilization operator being λ h the jump of λ h across the internal faces of G h . Then, we use the result of [7], Section III to show that which combined with (4.17) shows that the second assumption of Lemma 4.8 holds true. The third step of the analysis consists of showing that (4.10) and (4.11) are satisfied. We introduce the following discrete norms recalling that h is the mesh size of T Ω h . We equip the space X h with the discrete norm where δ is an arbitrarily small parameter. Then the following result holds true.
Proof. Concerning the coercivity property (4.10), we show that ∀[u h , u h ], there exists ξ h ∈ Q h such that, Using the definitions of π H and ξ h ([u h , u h ]) previously presented and recalling that ξ h ∈ Q H ⊂ Q h , we obtain Now, we seek for an upper bound of the second and third (negative) terms of the last inequality. For the second term, we apply the additional assumption that the operators T Λ and ∂ s commute. This is true if the cross section D does not depend on the arclength s. Then, we use the Poincaré inequality (4.19), the average inequality (4.13) and the trace inequality (4.12) to show that, For the third term, the following upper bound holds true, Combining the last three inequalities, reminding that we obtain the desired inequality. Concerning inequality (4.11), the proof is analogous to the one in [7].

A benchmark problem with analytical solution.
Let Ω = where for formulation (2.4) the mix-dimensional coupling constraint reads while for (2.2) we set In (5.1)-(5.3) the right-hand sides shall be defined as f = 8π 2 sin(2πx) sin(2πy), g = π 2 sin (πz) , u b = sin(2πx) sin(2πy), The exact solution of (5.1), regardless of the coupling constraint, is given by Let us notice that u satisfies homogeneous Dirichlet conditions at the boundary of Λ. Moreover, the solution (5.4) satisfies on Γ the relation (5.5) ∇u · n ⊕ = d z u n ⊕,z = 0, with n ⊕,z the z−component of the normal unit vector to Γ. We prove that the solution of (5.1) is equivalent to the one of (2.2). Precisely, we prove that (5.4) is the solution of (2.2). Using the integration by part formula and homogeneous boundary conditions on Ω and Λ, from (2.2) we have Since λ = 0 and the first of (5.4) satisfies (5.1a) and the second satisfies (5.1b), we have that Thus (5.4) satisfy equations (2.2a), (2.2b). The fact that the solution satisfy (2.2c) follows from (5.3). We can prove in a similar way that (5.4), with λ = 0 satisfy (2.4). Note in particular that q 1 is such that T Γ u − E Γ u = q 1 on Γ. T

Numerical experiments.
T Ω h conforming to Γ. Using the benchmark problem (5.1) we now investigate convergence properties of the two formulations. To this end we consider a uniform mesh of T Ω h of Ω consisting of tetrahedra with diameter h. Further, the discretization shall be geometrically conforming to both Λ and Γ such that the meshes T Γ h , T Λ h are made up of facets and edges of T Ω h respectively, cf. Figure  5.1 for illustration.
Considering inf-sup stable discretization in terms of continuous linear Lagrange (P 1 ) elements (for all the spaces), Table 5.1 lists the errors of formulations (2.4) and (2.2) on the benchmark problem. It can be seen that the error in u and u in H 1 norm converges linearly (as can be expected due to P 1 element discretization). Moreover, the error of the Lagrange multiplier approximation in H −1/2 norm decreases quadratically. In the light of P 1 discretization this rate appears superconvergent. We speculate that the result is due to the fact that the exact solution is particularly simple, λ = λ = 0. We remark that for u and u the error is interpolated into the finite element space of piecewise quadratic discontinous functions. For (2.2) we evaluate the fractional norm and interpolate the error using piecewise continuous cubic functions. For the sake of comparison with non-conforming formulation of (2.2) from  Table 5.1 also lists the error of the Lagrange multiplier in the L 2 norm. Here, quadratic convergence is observed for (2.2). For (2.4) the rate is between 1.5 and 2.
We plot the numerical solution of problem (2.4) and (2.2) in Figure 5.2.

Numerical experiments.
T Ω h non-conforming to Γ. Using benchmark problem (5.1) we consider (2.2) in the setting of §4.2. To this end we let T Ω h be a uniform mesh of Ω such that no cell T Ω h has any edge lying on Λ. Further we let h = h/3 in T Λ h , cf. Figure 5.1. Using discretization in terms of P 1 -P 1 -P 0 element Table 5.1 lists the error of the stabilized formulation of (2.2). A linear convergence in the H 1 norm can be observed in the error of u and u . We remark that the norms were computed as in §5.1. For simplicity the convergence of the multiplier is measured in the L 2 norm rather then the H −1/2 (Γ) norm used in the analysis. Then, convergence exceeding order 1.5 can be observed, however, the rates are rather unstable. Tables 5.1 one can observe that all the formulations yield practically identically accurate approximations of u. Further, compared to the conforming case, the stabilized formulation (2.2) results in a greater accuracy of u h as the underlying mesh T Λ h is here finer. Due to the different definitions in the three formulations, comparision of the Lagrange multiplier convergence is not straightforward. We therefore limit ourselves to a comment that in the L 2 norm all the formulations yield faster than linear convergence. In order to discuss solution cost of the formulations we consider the resulting preconditioned linear systems. In particular, we shall compare spectral condition numbers and the time to convergence of the preconditioned minimal residual (MinRes) solver with the with stopping criterion requiring the relative preconditioned residual norm to be less than 10 −8 . We remark that we shall ignore the setup cost of the preconditioner. Following operator preconditioning technique [23] we propose as preconditioners for (2.4) and (2.2) in the conforming case the (approximate) Riesz mapping with respect to the inner products of the spaces in which the two formulations were proved to be well posed. In particular, the preconditioner for the Lagrange multiplier relies on (the inverse of) the fractional Laplacian −∆ −1/2 on Γ for (2.4) and Λ for (2.2). A detailed analysis of the preconditioners will be presented in a separate work. We remark that in both cases the fractional Laplacian was here realized by spectral decomposition [21]. For the unfitted stabilized formulation (2.2) the Lagrange multiplier preconditioner uses a Riesz map with respect to the inner product due to L 2 (G h ) and the stabilization (4.20), i.e.  Table 5.2 Cost comparison of the formulations across refinement levels l. Number of Krylov iterations (preconditioned conjugate gradient for (5.1a), MinRes otherwise) and the condition number of the preconditioned problem is denoted by # and κ respectively. Time till convergence of the iterative solver (excluding the setup) is shown as T .

Comparison. In
This simple choice does not yield bounded iterations. However, establishing a robust preconditioner in this case is beyond the scope of the paper and shall be pursued in the future works. In Table 5.2 we compare solution time, number of iterations and condition numbers of the (linear systems due to the) three formulations. Let us first note that the proposed preconditioners for (2.4) and (2.2) in the conforming case seem robust with respect to discretization parameter as the iteration counts and condition numbers are bounded in h. We then see that the solution time for (2.4) is about 2 times longer compared to (2.2) which is about 4 times more expensive than the solution of the Poisson problem (5.1a). This is in addition to the higher setup costs of the preconditioner, which in our implementation involve solving an eigenvalue problem for the fractional Laplacian. Therefore it is advantageous to keep the multiplier space as small as possible. We remark that the missing results for (2.4) in Table 5.2 are due to the memory limitations encountered when solving the eigenvalue problem for the Laplacian, which for finest mesh involves cca 32 thousand eigenvalues, cf. Appendix C. Due to the missing proper preconditioner for the Lagrange multiplier block the number of iterations in the third, unfitted formulation can be seen to approximately double on refinement.
Appendix A. Derivation of the model. This section provides a rigorous derivation of 3D-1D-1D problem (2.2) and 3D-1D-2D problem (2.4). The steps are similar to the derivation presented in [22], however, here the coupling conditions are different, giving rise to coupled problems featuring Lagrange multipliers. Precisely, the starting point is the problem arising from Dirichlet-Neumann conditions. Find u ⊕ , u s.t.: The coupling constraints defined on Γ involve essential or strong conditions. Such conditions will be enforced weakly by using the method of Lagrange multipliers [2]. Then, the variational formulation of problem (A.1) is to find u ⊕ ∈ H 1 ∂Ω (Ω ⊕ ), u ∈ H 1 ∂Ω \Γ (Ω ), λ ∈ H − 1 2 (Γ) s.t.
where λ is the Lagrange multiplier and it is equivalent to ∇u · n .
Model reduction of the problem on Ω . We apply the averaging technique to equation (A.1b). In particular, we consider an arbitrary portion P of the cylinder Ω , with lateral surface Γ P and bounded by two perpendicular sections to Λ, namely D(s 1 ), D(s 2 ) with s 1 < s 2 . We have, From the combination of all the above terms with the right hand side, we obtain that the solution u of (A.1b) satisfies, Let us now formulate the modelling assumption that allows us to reduce equation (A.5) to a solvable one-dimensional (1D) model. More precisely, we assume that the function u has a uniform profile on each cross section D(s), namely u (r, s, t) = u (s). Therefore, observing that u = u = u , and that ∂ s u = ∂ s u = d s u , problem (A.5) turns out to find u ∈ H 1 (Λ) such that (A.6) (d s u , d s v ) Λ,|D| +(u , v ) Λ,|D| −(λ, v ) Λ,|∂D| = (g, v ) Λ,|D| ∀v ∈ H 1 (Λ) .
Topological model reduction of the problem on Ω ⊕ . We focus here on the subproblem of (A.1a) related to Ω ⊕ . We multiply both sides of (A.1a) by a test function v ∈ H 1 0 (Ω) and integrate on Ω ⊕ . Integrating by parts and using boundary and interface conditions, we obtain Then, we make the following modelling assumption: we identify the domain Ω ⊕ with the entire Ω, and we correspondingly omit the subscript ⊕ to the functions defined on Ω ⊕ , namely Ω⊕ u ⊕ dω Ω u dω .