Coupled partial differential equations (PDEs) defined on domains with different dimensionality are usually called mixed-dimensional PDEs. We address mixed-dimensional PDEs on three-dimensional (3D) and one-dimensional (1D) domains, which gives rise to a 3D-1D coupled problem. Such a 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

more »

... ns, which are 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 of determining 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. 559 the coupling constraint (1.1c); see section SM1 of the supplementary material for a precise definition. 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 km [27, 28] . Similarly, in applications involving blood flow and oxygen transport of microcirculation, the capillary radius is a few microns, while simulations are often performed on mm to cm scale, with thousands of vessels [3, 14, 17, 32] . Finally, in neuroscience applications a neuron has a width of a few microns, while its length is much longer. For example, an axon of a motor neuron may be as long as a meter. Hence, at least four orders of magnitude in difference in the transverse and longitudinal directions is common in geophysics, biomechanics, and neuroscience. Meshes dictated by resolving the transverse length scale in three dimensions would then possibly lead to the order of 10 12 degrees of freedom. Even if adaptive and strongly anisotropic meshes were allowed, the computations would quickly become demanding if many slender structures and their interactions were considered. From a mathematical standpoint, the challenge involved in problem (1.1) is that neither \scrT \Lambda nor \delta \Lambda is well defined. That is, without extra regularity, solutions of elliptic PDEs only have well defined traces of codimension one. Here, \scrT \Lambda is of codimension two, mapping functions defined on a domain in three dimensions 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 works by D'Angelo [10, 11] and D'Angelo and Quarteroni [12] have remedied the well-posedness by weakening the solution concept. This approach naturally leads to nonsymmetric formulations. An alternative approach is to decompose the solution into smooth and nonsmooth components, where the nonsmooth component may be represented in terms of Green's functions, and then consider the well-posedness of the smooth component [16] . 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 that has been addressed in [4, 18, 20] , where the optimal a priori error estimates for the finite element approximation are derived. Furthermore, the interplay among the mathematical structure of the problem, its solvers, and preconditioners for its discretization was studied in detail in [22] for the solution of 1D differential equations embedded in two dimensions and more recently was extended to the 3D-1D case in [21] . Stemming from this literature, our work adopts and analyzes a different approach closely related to [19, 23] . That is, we exploit the fact that \Lambda is not strictly a 1D curve but rather a very thin 3D structure with a cross-sectional area far below 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 of this work is that we address essential coupling conditions, namely Dirichlet--Neumann conditions; see, in particular, problem (SM1.1) in the supplementary material. In previous works, for example, [12, 19, 23] , natural coupling conditions of Robin--Robin type were analyzed. Dirichlet-type coupling conditions pose additional difficulties as the conditions are not a natural part of the weak formulation of the problem. As shown in section SM1 of the supplementary material, 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 on the analysis and approximation of the proposed approach, we stress that it aims to build the mathematical foundations for tackling various applications involving 3D-1D mixed-dimensional PDEs, such as fluid-structure interaction of slender bodies [26] , microcirculation and lymphatics [29, 33] , subsurface flow models with wells [8] , and the electrical activity of neurons. Preliminaries. Let the domain \Omega \subset \BbbR 3 be an open, connected, and convex set that can be subdivided into two parts \Omega \ominus and \Omega \oplus := \Omega \setminu \Omega \ominus . Let \Omega \ominus be a generalized cylinder (cf. [15] ), that is, the swept volume of a two-dimensional set, \partial\scrD , moved along a curve, \Lambda , in the three-dimensional domain, \Omega ; see Figure 2 .1 for an illustration. More precisely, the curve \Lambda = \{ \bfitlam (s), s \in (0, S)\} , where \bfitlam (s) = [\xi (s), \tau (s), \zeta (s)], and s \in (0, S) is a \scrC 2 -regular curve in the three-dimensional domain \Omega . For simplicity, let us assume that \| \bfitlam \prime (s)\| = 1 such that the arclength and the coordinate s coincide. Further, let \scrD (s) = [x(r, t), y(r, t)] : (0, R(s)) \times (0, T (s)) \rightar \BbbR 2 be a parametrization of the cross-section with R(s) \geq R 0 > 0 being R 0 , that is, the minimum cross-sectional radius of the generalized cylinder, and let \Gamma be the lateral surface of \Omega \ominus , i.e., \Gamma = \{ \partial\scrD (s) | s \in \Lambda \} , while the upper and lower faces of \Omega \ominus belong to \partial\Omega . We assume that \Omega \ominus crosses \Omega from side to side. Finally, | \cdot | denotes the Lebesgue measure of a set, e.g., | \scrD (s)| is the cross-sectional area of the cylinder. In general, | \scrD (s)| must be strictly positive and bounded. According to the geometrical setting, we will denote by v, v \oplus , v \ominus , v \odot functions defined on \Omega , \Omega \oplus , \Omega \ominus , \Lambda , respectively. Let D be a generic regular bounded domain in \BbbR 3 , and let X be a Hilbert space defined on D. Then (\cdot , \cdot ) X and \| \cdot \| X denote the inner product and norm of X, respectively. The duality pairing between the X and its dual X \ast is denoted by \langle \cdot , \cdot \rangle . Let (\cdot , \cdot ) L 2 (D) , (\cdot , \cdot ) D , or simply (\cdot , \cdot ) 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 \| \cdot \| H q (D) , and the seminorm is | \cdot | 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 \Sigma be a Lipschitz codimension one subset of D. We denote by \scrT \Sigma : H q (D) \rightar H q -1 2 (\Sigma ) the trace operator from D to \Sigma . The space of functions in H 1 2 00 (\Lambda ). 3. Saddle-point problem analysis. Let a : X \times X \rightar \BbbR and b : X \times Q \rightar \BbbR be bilinear forms. Let us consider a general saddle-point problem of the following form: Downloaded 03/10/21 to 129.240.222.175. Redistribution subject to SIAM license or copyright; see https://epubs.siam.org/page/terms