3D-1D coupling on non conforming meshes via three-field optimization based domain decomposition
33D-1D coupling on non conforming meshes via three-fieldoptimization based domain decomposition
Stefano Berrone a, ∗ , Denise Grappein a , Stefano Scialò a a Dipartimento di Scienze Matematiche, Politecnico di Torino, Corso Duca degli Abruzzi 24,10129 Torino, Italy. Member of INdAM research group GNCS.
Abstract
A new numerical approach is proposed for the simulation of coupled three-dimensional and one-dimensional elliptic equations (3D-1D coupling) arisingfrom dimensionality reduction of 3D-3D problems with thin inclusions. Themethod is based on a well posed mathematical formulation and results in anumerical scheme with high robustness and flexibility in handling geometricalcomplexities. This is achieved by means of a three-field approach to split the1D problems from the bulk 3D problem, and then resorting to the minimizationof a properly designed functional to impose matching conditions at the inter-faces. Thanks to the structure of the functional, the method allows the use ofindependent meshes for the various subdomains.
Keywords:
1. Introduction
This work presents a new numerical approach to manage the coupling ofthree-dimensional and one-dimensional elliptic equations (3D-1D coupling). Thiskind of problems emerges, for example, in the numerical treatment of domainswith small tubular inclusions: in these cases, indeed, it might result computa-tionally convenient to approximate the small inclusions by one-dimensional (1D)manifolds, in order to avoid the building of a three-dimensional grid within theinclusion. Clearly this topological reduction can be a viable approach only ifone-dimensional modeling assumptions can be applied to the problem at hand.Examples of application are: capillary networks exchanging flux with the sur-rounding tissue [1, 2], the interaction of tree roots with the soil [3, 4], a systemof wells for fluid in geological applications [5, 6, 7], or the modeling of fiber-reinforced materials [8, 9]. ∗ Corresponding author
Email addresses: [email protected] (Stefano Berrone ), [email protected] (Denise Grappein), [email protected] (Stefano Scialò)
Preprint submitted to Elsevier February 15, 2021 a r X i v : . [ m a t h . NA ] F e b he definition of coupling conditions between a three-dimensional (3D) anda 1D problem is not straightforward, as no bounded trace operator is definedin standard functional spaces on manifolds with a dimensionality gap higherthan one. This problem was recently studied in [10, 11], where suitable weighedSobolev spaces were introduced and a bounded trace operator from the 3D spaceto the 1D space was defined, thus allowing to formulate a well posed coupledproblem, also resorting to the results in [12]. In [13] regularizing techniques areproposed for singular terms. Three dimensional problems with singular sourcesdefined on lines are also studied in [14], where the nature of the irregularity isanalyzed and a method based on the splitting of the solution in a low regularitypart plus a regular correction is proposed. Problems with a source term onmanifolds with high dimensionality gaps are also studied in [15], where a liftingtechnique of the irregular datum is used to reduce the dimensionality gap. In[16] a 3D-1D coupled approach is derived starting from the fully 3D-3D coupledproblem and applying a topological model reduction through the definition ofproper averaging operators.In the present work a new numerical approach to this problem is proposed,starting from modeling assumptions similar to those proposed in [16]. A re-duced 3D-1D model approximating the original equi-dimensional 3D-3D prob-lem is here obtained by introducing proper assumptions on the solution insidethe small inclusions and by defining suitable subspaces of the Sobolev spacestypically employed in the variational formulation of partial differential equa-tions. Thanks to this, the problem is treated, in practice, as a 3D-1D reducedproblem, but it can still be written as an equi-dimensional problem, thus skip-ping the difficulties related to the 3D-1D coupling. The problems in the bulk3D domain and in the small inclusions are splitted resorting to a three-fieldbased domain decomposition method, originally formulated in [17] and alreadyapplied for domain decomposition in networks of fractures [18]. Suitable match-ing conditions are then enforced at the interface to recover the solution on thewhole domain. Here pressure continuity and flux conservation constraints areassumed at the interfaces. The advantages of the three-field based approachwith respect to standard domain decomposition methods lie in the possibilityof defining stable locally conservative numerical schemes on non conformingmeshes [18]. Matching conditions at the interfaces are enforced by means of aPDE constrained numerical scheme, already proposed for simulation of the flowin poro-fractured media [19, 20, 21] and here proposed for 3D-1D coupling. Themethod is based on the minimization of a cost functional expressing the error inthe fulfillment of interface conditions, constrained by constitutive laws on thevarious subdomains. The advantages of such an approach lie in the possibilityof enforcing continuity and flux balance at the interfaces using completely inde-pendent meshes on the sub-domains. This provides to the scheme an extremeflexibility and robustness to geometrical complexities, which are critical featuresfor the applicability to real problems, where the nearly one dimensional inclu-sions might form complex networks. Exploiting the properties of the functionalspaces chosen for the solution, the functional can be reduced to the centrelines ofthe tubular inclusions and used to control the continuity of the solution, whereas2ux conservation is strongly enforced thanks to the three-field formulation.The manuscript is organized as follows: in Section 2 notation is introducedand the strong formulation of the 3D-3D problem is presented, along with thehypotheses allowing its reduction to a 3D-1D problem. In Section 3 the weakformulation of the 3D-1D coupled problem is worked out, while in Section 4 theproblem is re-written into a PDE-constrained optimization formulation. Thecorresponding discrete approach is discussed in Section 6. Finally, in Section 7,some numerical examples are described.
2. Notation and problem formulation
Let us here briefly recall the basic formulation of the problem of interest ina simplified setting, the ideas here proposed being easily extendable to moregeneral cases. We refer to [16] for a broader presentation of the problem. Letus consider a convex domain Ω ∈ R in which a generalized cylinder Σ ∈ R isembedded. The centreline of this cylinder is denoted by Λ = { λ ( s ) , s ∈ (0 , S ) } ,where λ ( s ) is here assumed, for simplicity, to be a rectilinear segment in thethree-dimensional space. The symbol Σ( s ) denotes the transversal section ofthe cylinder at s ∈ [0 , S ] . We assume that each section, whose boundary isdenoted by Γ( s ) , has an elliptic shape, denoting by R the maximum axes lengthof the ellipses as s ranges in the interval [0 , S ] . The lateral surface of Σ is Γ = { Γ( s ) , s ∈ [0 , S ] } , while Σ = Σ(0) and Σ S = Σ( S ) are the two extremesections. The portion of the domain that does not include the cylinder is denotedby D = Ω \ Σ , with boundary ∂D = ∂ Ω ∪ { Γ ∪ Σ ∪ Σ S } , where ∂ Ω is theboundary of Ω . We refer to ∂D e as the external boundary of D , coinciding with ∂ Ω when the extreme sections of Σ are inside Ω . In case Σ and Σ S lie on theboundary ∂ Ω we define ∂D e = ∂ Ω \ { Σ ∪ Σ S } .Let us consider, in Ω , a diffusion problem, with unknown pressures u in D and ˜ u in Σ : D : −∇ · ( K ∇ u ) = f in D (1) u = 0 on ∂D e (2) u | Γ = ψ on Γ (3) K ∇ u · n = φ on Γ (4) Σ : −∇ · ( ˜ K ∇ ˜ u ) = g in Σ (5) ˜ u = 0 on Σ ∪ Σ S (6) ˜ u | Γ = ψ on Γ (7) ˜ K ∇ ˜ u · ˜ n = − φ on Γ (8)Vectors n and ˜ n are the outward pointing unit normal vectors to Γ , respectivelyfor D and Σ , such that ˜ n = − n , K and ˜ K are uniformly positive definitetensors in D and Σ , respectively, and f and g are source terms. The symbol ψ denotes the unknown unique value of the pressure on the interface Γ while φ isthe unknown flux through Γ , entering in D . Homogeneous Dirichlet boundaryconditions are considered on ∂D e and on Σ and Σ S , assumed to be part of ∂ Ω ,the extension to more general cases being straightforward. Equations (3)-(4)and (7)-(8) enforce pressure continuity and flux conservation constraints on thelateral surface Γ of the cylinder, thus allowing us to couple the two problems.3e remark that different coupling conditions could be considered, as the onesproposed in [16].If the cross-section-size R of the cylinder Σ becomes much smaller than thecharacteristic dimension of Ω , a model can be introduced in order to reduce thecomputational cost of simulations. The key point is that, as R (cid:28) diam (Ω) , itis possible to assume that the variations of ˜ u on the transversal sections of thecylinder are negligible, i.e., in cylindrical coordinates, ˜ u ( r, θ, s ) = ˆ u ( s ) ∀ r ∈ [0 , R ] , ∀ θ ∈ [0 , π ) . (9)This allows us to simplify problem (5)-(6) reducing it to a 1D problem definedon the cylinder’s centerline Λ as − dds (cid:32) ˜ K | Σ( s ) | d ˆ uds (cid:33) = ˜ g for s ∈ (0 , S ) (10) ˆ u (0) = ˆ u ( S ) = 0 , (11)where the new forcing term ˜ g now accounts for the original volumetric source g and for the incoming flux from the boundary Γ of the equi-dimensional problem.As the domain of problem (5)-(6) is reduced to a 1D segment, the domain ofproblem (1)-(2) is extended to Ω \ Λ , thus obtaining a coupled 3D-1D problem.Details of this geometrical reduction are provided in the next session. Theclear advantage of the reduced problem is that solving a problem defined ona segment instead of a problem defined on a small cylinder is computationallymuch cheaper. Nevertheless it is not straightforward to write coupling conditionsbetween (10)-(11) and (1)-(2) that are analogous to (3)-(4) and (7)-(8), as thereis no bounded trace operator from H (Ω) to L (Λ) , being the dimensionalitygap between the interested manifolds higher than one.In the next section an original formulation of the coupled 3D-1D problemis derived, starting from a variational formulation of the fully dimensional 3D-3D problem (1)-(2) and (5)-(6), and defining the proper functional spaces andoperators required to reduce this formulation to a well posed 3D-1D coupling.
3. Variational formulation
In this section we will adopt the following notation: H ( D ) = (cid:8) v ∈ H ( D ) : v | ∂De = 0 (cid:9) H (Σ) = (cid:110) v ∈ H (Σ) : v | Σ0 = v | Σ S = 0 (cid:111) H (Λ) = (cid:8) v ∈ H (Λ) : v (0) = v ( S ) = 0 (cid:9) and we will suppose that ∂D e = ∂ Ω \ { Σ ∪ Σ S } , i.e. the extreme sections ofthe cylinder lie on ∂ Ω . Let us define a trace operator γ Γ : H ( D ) ∪ H (Σ) → H (Γ) , such that γ Γ v = v | Γ ∀ v ∈ H ( D ) ∪ H (Σ) (12)4nd two extension operators E Σ : H (Λ) → H (Σ) and E Γ : H (Λ) → H (Γ) defined such that, given ˆ v ∈ H (Λ) , E Σ (ˆ v ) is the extension of the point-wise value ˆ v ( s ) , s ∈ [0 , S ] , to the cross section Σ( s ) of the cylinder and E Γ (ˆ v ) the extensionto the boundary Γ( s ) of the cross section. Let us observe that E Γ = γ Γ ◦ E Σ .Setting ˆ V = H (Λ) , let us further consider the spaces: (cid:101) V = { v ∈ H (Σ) : v = E Σ ˆ v, ˆ v ∈ ˆ V } , H Γ = { v ∈ H (Γ) : v = E Γ ˆ v, ˆ v ∈ ˆ V } ,V D = (cid:8) v ∈ H ( D ) : γ Γ v ∈ H Γ (cid:9) . such that (cid:101) V ⊂ H (Σ) contains functions that are extensions to the whole Σ of functions in ˆ V , H Γ ⊂ H (Γ) contains functions that are extensions to Γ offunctions in ˆ V , and V D ⊂ H ( D ) only contains functions whose trace on Γ is afunction of H Γ . The variational problem arising from the coupling of (1)-(2) and(5)-(6) through the continuity constraint can be written as: find ( u, ˜ u ) ∈ V D × (cid:101) V such that ( K ∇ u, ∇ v ) V D + ( ˜ K ∇ ˜ u, ∇ ˜ v ) (cid:101) V = ( f, v ) V D + ( g, ˜ v ) (cid:101) V ∀ ( v, ˜ v ) ∈ V D × (cid:101) V (13) (cid:104) γ Γ u − γ Γ ˜ u, η (cid:105) H Γ , H Γ (cid:48) = 0 , ∀ η ∈ H Γ (cid:48) (14) Remark 1.
Let us consider the space V = (cid:110) ( v, ˜ v ) ∈ V D × (cid:101) V : γ Γ v = γ Γ ˜ v (cid:111) .Then, problem (13)-(14) is equivalent to: find ( u, ˜ u ) ∈ V such that ( K ∇ u, ∇ v ) V D + ( ˜ K ∇ ˜ u, ∇ ˜ v ) (cid:101) V = ( f, v ) V D + ( g, ˜ v ) (cid:101) V ∀ ( v, ˜ v ) ∈ V The well-posedness of the problem easily follows from Lax-Milgram theorem,considering || · || V = || · || H ( D ) + || · || H (Σ) Equation (13) can be split into two coupled equations introducing the un-known flux φ through Γ , obtaining thus ( K ∇ u, ∇ v ) V D − (cid:104) φ, γ Γ v (cid:105) H Γ (cid:48) , H Γ = ( f, v ) V D ∀ v ∈ V D , φ ∈ H Γ (cid:48) (15) ( ˜ K ∇ ˜ u, ∇ ˜ v ) (cid:101) V + (cid:104) φ, γ Γ ˜ v (cid:105) H Γ (cid:48) , H Γ = ( g, ˜ v ) (cid:101) V ∀ ˜ v ∈ (cid:101) V , φ ∈ H Γ (cid:48) (16)Moreover, the continuity condition (14) can be rewritten introducing a newvariable ψ ∈ H Γ as: (cid:104) γ Γ u − ψ, η (cid:105) H Γ , H Γ (cid:48) = 0 ∀ η ∈ H Γ (cid:48) , ψ ∈ H Γ (17) (cid:104) γ Γ ˜ u − ψ, η (cid:105) H Γ , H Γ (cid:48) = 0 ∀ η ∈ H Γ (cid:48) , ψ ∈ H Γ . (18)The set of equations (15), (16), (17) and (18) represent an application of thethree-field formulation presented in [17] and similarly applied in [18].5hanks to the assumptions on the introduced functional spaces, this problemcan be reduced to a 3D-1D coupled problem without encountering the aforemen-tioned issues in the definition of a trace operator. In fact we only need to usethe trace operator γ Γ ( · ) , which is well-posed as defined from a three-dimensionalmanifold to a two dimensional one. Let us observe that (cid:104) φ, γ Γ v (cid:105) H Γ (cid:48) , H Γ = (cid:90) Γ φ γ Γ v d Γ = (cid:90) S (cid:16) (cid:90) Γ( s ) φ γ Γ v dl (cid:17) ds ∀ v ∈ V D and let us denote by φ ( s ) the mean value of φ on the border Γ( s ) of eachsection. As v ∈ V D we know that γ Γ v ∈ H Γ , i.e. ∃ ˇ v ∈ ˆ V : γ Γ v = E Γ ˇ v . Thus (cid:82) Γ( s ) γ Γ v dl = | Γ( s ) | ˇ v ( s ) and (cid:90) S (cid:16) (cid:90) Γ( s ) φ γ Γ v dl (cid:17) ds = (cid:90) S | Γ( s ) | φ ( s )ˇ v ( s ) ds = (cid:10) | Γ | φ, ˇ v (cid:11) ˆ V (cid:48) , ˆ V , where | Γ( s ) | is the section perimeter size at s ∈ [0 , S ] . The same holds if weconsider (cid:104) φ, γ Γ ˜ v (cid:105) H Γ (cid:48) , H Γ . As ˜ v ∈ (cid:101) V we know that ∃ ˆ v ∈ ˆ V : ˜ v = E Σ ˆ v andconsequently γ Γ ˜ v = γ Γ E Σ ˜ v = E Γ ˆ v , so that (cid:104) φ, γ Γ ˜ v (cid:105) H Γ (cid:48) , H Γ = (cid:10) | Γ | φ, ˆ v ( s ) (cid:11) ˆ V (cid:48) , ˆ V . Similarly, ∀ η ∈ H Γ (cid:48) and ∀ ρ : ρ = E Γ ˆ ρ with ˆ ρ ∈ ˆ V , we can write (cid:104) ρ, η (cid:105) H Γ (cid:48) , H Γ = (cid:90) S (cid:16) (cid:90) Γ( s ) ρηdl (cid:17) ds = (cid:90) S | Γ( s ) | ˆ ρ ( s ) η ( s ) ds = (cid:104) ˆ ρ, | Γ | η (cid:105) ˆ V (cid:48) , ˆ V (19)where we have used (cid:82) Γ( s ) ρ dl = | Γ( s ) | ˆ ρ and η ( s ) is the mean value of η on theborder Γ( s ) of each section. Exploiting (19) we can rewrite conditions (17) and(18) as (cid:104) γ Γ u − ψ, η (cid:105) H Γ , H Γ (cid:48) = (cid:68) | Γ | (ˇ u − ˆ ψ ) , η (cid:69) ˆ V , ˆ V (cid:48) = 0 (cid:104) γ Γ ˜ u − ψ, η (cid:105) H Γ , H Γ (cid:48) = (cid:68) | Γ | (ˆ u − ˆ ψ ) , η (cid:69) ˆ V , ˆ V (cid:48) = 0 where ˇ u, ˆ ψ ∈ ˆ V are such that γ Γ u = E Γ ˇ u , ψ = E Γ ˆ ψ and γ Γ ˜ u = γ Γ E Σ ˜ u = E Γ ˆ u , as ˜ u ∈ (cid:101) V . Finally let us observe that ( ˜ K ∇ ˜ u, ∇ ˜ v ) (cid:101) V = (cid:90) Σ ˜ K ∇ ˜ u ∇ ˜ v dσ = (cid:90) S ˜ K | Σ( s ) | d ˆ uds d ˆ vds ds where ˆ u, ˆ v ∈ ˆ V are such that ˜ u = E Σ ˆ u , ˜ v = E Σ ˆ v and | Σ( s ) | is the sectionarea at s ∈ [0 , S ] . If we now extend the space V D from D to the whole region Ω , denoting this extended space by V , problem (15)-(18) can be reduced to a6D-1D coupled problem as: Find ( u, ˆ u ) ∈ V × ˆ V , φ ∈ ˆ V (cid:48) and ˆ ψ ∈ ˆ V such that: ( K ∇ u, ∇ v ) V − (cid:10) | Γ | φ, ˇ v (cid:11) ˆ V (cid:48) , ˆ V = ( f, v ) V ∀ v ∈ V, ˇ v ∈ ˆ V : γ Γ v = E Γ ˇ v (20) (cid:16) ˜ K | Σ | d ˆ uds , d ˆ vds (cid:17) ˆ V + (cid:10) | Γ | φ, ˆ v (cid:11) ˆ V (cid:48) , ˆ V = ( | Σ | g, ˆ v ) ˆ V ∀ ˆ v ∈ ˆ V (21) (cid:68) | Γ | (ˇ u − ˆ ψ ) , η (cid:69) ˆ V (cid:48) , ˆ V = 0 γ Γ u = E Γ ˇ u, ∀ η ∈ ˆ V (cid:48) (22) (cid:68) | Γ | (ˆ u − ˆ ψ ) , η (cid:69) ˆ V (cid:48) , ˆ V = 0 ∀ η ∈ ˆ V (cid:48) (23)with g ( s ) = | Σ( s ) | (cid:82) Σ( s ) g dσ , being g sufficiently regular.
4. PDE-constrained optimization problem
The fulfillment of conditions (22) and (23) can be obtained through theminimization of a cost functional. Since we want to formulate independentproblems on the various sub-domains, in order to guarantee the well posedness ofeach problem independently from the imposed boundary conditions, we modifyequations (20)-(21) as follows: ( K ∇ u, ∇ v ) V + α ( | Γ | ˇ u, ˇ v ) ˆ V − (cid:10) | Γ | φ, ˇ v (cid:11) ˆ V (cid:48) , ˆ V = ( f, v ) V + α ( | Γ | ˆ ψ, ˇ v ) ˆ V (24) ∀ v ∈ V, ˇ v ∈ ˆ V : γ Γ v = E Γ ˇ v, (cid:16) ˜ K | Σ | d ˆ uds, d ˆ vds (cid:17) ˆ V + ˆ α ( | Γ | ˆ u, ˆ v ) ˆ V + (cid:10) | Γ | φ, ˆ v (cid:11) ˆ V (cid:48) , ˆ V = ( | Σ | g, ˆ v ) ˆ V + ˆ α ( | Γ | ˆ ψ, ˆ v ) ˆ V (25) ∀ ˆ v ∈ ˆ V . with α, ˆ α > being arbitrary parameters. Let us now define the followingfunctional J ( φ, ˆ ψ ) = 12 (cid:16) || γ Γ u ( φ, ˆ ψ ) − ψ || H Γ + || γ Γ ˜ u ( φ, ˆ ψ ) − ψ || H Γ (cid:17) = 12 (cid:16) || γ Γ u ( φ, ˆ ψ ) − E Γ ˆ ψ || H Γ + || γ Γ E Σ ˆ u ( φ, ˆ ψ ) − E Γ ˆ ψ || H Γ (cid:17) (26)to be minimized constrained by (24) and (25). In order to rewrite the PDE-constrained optimization problem in a compact form, we consider the linearoperators A : V → V (cid:48) , (cid:98) A : ˆ V → ˆ V (cid:48) , B : ˆ V (cid:48) → V (cid:48) , (cid:98) B : ˆ V (cid:48) → ˆ V (cid:48) , C : ˆ V → V (cid:48) and (cid:98) C : ˆ V → ˆ V (cid:48) such that: (cid:104) Au, v (cid:105) V (cid:48) ,V = ( K ∇ u, ∇ v ) V + α ( | Γ | ˇ u, ˇ v ) ˆ V v ∈ V, ˇ v ∈ ˆ V : γ Γ v = E Γ ˇ v (27) (cid:68) (cid:98) A ˆ u, ˆ v (cid:69) ˆ V (cid:48) , ˆ V = (cid:16) ˜ K | Σ | d ˆ uds, d ˆ vds (cid:17) ˆ V + ˆ α ( | Γ | ˆ u, ˆ v ) ˆ V ˆ v ∈ ˆ V (28)7 Bφ, v (cid:11) V (cid:48) ,V = (cid:10) | Γ | φ, ˇ v (cid:11) ˆ V (cid:48) , ˆ V v ∈ V, ˇ v ∈ ˆ V : γ Γ v = E Γ ˇ v (29) (cid:68) (cid:98) Bφ, ˆ v (cid:69) ˆ V (cid:48) , ˆ V = (cid:10) | Γ | φ, ˆ v (cid:11) ˆ V (cid:48) , ˆ V ˆ v ∈ ˆ V (30) (cid:68) C ˆ ψ, v (cid:69) V (cid:48) ,V = α ( | Γ | ˆ ψ, ˇ v ) ˆ V v ∈ V, ˇ v ∈ ˆ V : γ Γ v = E Γ ˇ v (31) (cid:68) (cid:98) C ˆ ψ, ˆ v (cid:69) ˆ V (cid:48) , ˆ V = ˆ α ( | Γ | ˆ ψ, ˆ v ) ˆ V ˆ v ∈ ˆ V . (32)The respective adjoints will be denoted as A ∗ : V → V (cid:48) , (cid:98) A ∗ : ˆ V → ˆ V (cid:48) , B ∗ : V → ˆ V , (cid:98) B ∗ : ˆ V → ˆ V , C ∗ : V → ˆ V (cid:48) , (cid:98) C ∗ : ˆ V → ˆ V (cid:48) . Let us further define F ∈ V (cid:48) s.t. F ( v ) = ( f, v ) V , v ∈ V (33) G ∈ ˆ V (cid:48) s.t. G (ˆ v ) = ( | Σ | g, ˆ v ) ˆ V , ˆ v ∈ ˆ V . (34)Equations (24)-(25) can thus be written as: Au − Bφ − C ˆ ψ = F (35) (cid:98) A ˆ u + (cid:98) Bφ − (cid:98) C ˆ ψ = G. (36)If we now consider the space V = V × ˆ V and we set W = ( u, ˆ u ) ∈ V and V = ( v, ˆ v ) ∈ V , we can introduce the following operators: A : V → V (cid:48) s.t. A ( W , V ) = A ( u, v ) + (cid:98) A (˜ u, ˜ v ) B : ˆ V (cid:48) → V (cid:48) s.t. B ( φ, V ) = B ( φ, v ) − (cid:98) B ( φ, ˆ v ) C : ˆ V → V (cid:48) s.t. C ( ˆ ψ, V ) = C ( ˆ ψ, v ) + (cid:98) C ( ˆ ψ, ˆ v ) and the PDE-constrained optimization problem can be written as min ( φ, ˆ ψ ) J ( φ, ˆ ψ ) subject to (37) AW − B φ − C ˆ ψ = F , (38)with F ∈ V (cid:48) s.t. F ( V ) = F ( v ) + G (ˆ v ) . Proposition 1.
Let us consider the trace operator γ Γ : V → H Γ and the exten-sion operators E Σ : ˆ V → ˜ V ⊂ V and E Γ = γ Γ ◦ E Σ : ˆ V → H Γ , whose respectiveadjoints are γ ∗ Γ : H Γ (cid:48) → V (cid:48) , E Σ ∗ : ˜ V (cid:48) → ˆ V (cid:48) and E Γ ∗ : H Γ (cid:48) → ˆ V (cid:48) and let Θ ˆ V : ˆ V → ˆ V (cid:48) and Θ H Γ : H Γ → H Γ (cid:48) be Riesz isomorphisms. Then the optimalcontrol ( φ, ˆ ψ ) that provides the solution to (37) - (38) is such that Θ ˆ V ( B ∗ p − (cid:98) B ∗ ˆ p ) = 0 (39) Θ − V ( C ∗ p + (cid:98) C ∗ ˆ p − E Γ ∗ Θ H Γ ( γ Γ u ( φ, ˆ ψ ) + E Γ ˆ u ( φ, ˆ ψ ) − E Γ ˆ ψ )) = 0 (40) where p ∈ V and ˆ p ∈ ˆ V are the solutions respectively to A ∗ p = γ ∗ Γ Θ H Γ ( γ Γ u ( φ, ˆ ψ ) − E Γ ˆ ψ ) (41) (cid:98) A ∗ ˆ p = E Γ ∗ Θ H Γ ( E Γ ˆ u ( φ, ˆ ψ ) − E Γ ˆ ψ ) (42)8 roof. Let us compute the Frechet derivatives of the functional with respect tothe control variables φ and ˆ ψ . To this end, we introduce increments δφ ∈ ˆ V (cid:48) of φ and δ ˆ ψ ∈ ˆ V of ˆ ψ and we recall that ∃ δψ ∈ H Γ : δψ = E Γ δ ˆ ψ . We have: ∂J∂φ (cid:0) φ + δφ, ˆ ψ (cid:1) = (cid:16) γ Γ u ( φ, ˆ ψ ) − ψ, γ Γ u ( δφ, (cid:17) H Γ + (cid:16) γ Γ ˜ u ( φ, ˆ ψ ) − ψ, γ Γ ˜ u ( δφ, (cid:17) H Γ == (cid:16) γ Γ u ( φ, ˆ ψ ) − E Γ ˆ ψ, γ Γ u ( δφ, (cid:17) H Γ + (cid:16) γ Γ E Σ ˆ u ( φ, ˆ ψ ) − E Γ ˆ ψ, γ Γ E Σ ˆ u ( δφ, (cid:17) H Γ == (cid:68) γ ∗ Γ Θ H Γ ( γ Γ u ( φ, ˆ ψ ) − E Γ ˆ ψ ) , u ( δφ, (cid:69) V (cid:48) ,V + (cid:68) E Γ ∗ Θ H Γ ( E Γ ˆ u ( φ, ˆ ψ ) − E Γ ˆ ψ ) , ˆ u ( δφ, (cid:69) ˆ V (cid:48) , ˆ V == (cid:10) A ∗ p, A − Bδφ (cid:11) V (cid:48) ,V − (cid:68) (cid:98) A ∗ ˆ p, (cid:98) A − (cid:98) Bδφ (cid:69) ˆ V (cid:48) , ˆ V == (cid:10) B ∗ p, δφ (cid:11) ˆ V , ˆ V (cid:48) − (cid:68) (cid:98) B ∗ ˆ p, δφ (cid:69) ˆ V , ˆ V (cid:48) = (cid:16) Θ ˆ V ( B ∗ p − (cid:98) B ∗ ˆ p ) , δφ (cid:17) ˆ V (cid:48) ; ∂J∂ ˆ ψ ( φ, ˆ ψ + δ ˆ ψ ) = (cid:16) γ Γ u ( φ, ˆ ψ ) − ψ, γ Γ u (0 , δ ˆ ψ ) − δψ (cid:17) H Γ + (cid:16) γ Γ ˜ u ( φ, ˆ ψ ) − ψ, γ Γ ˜ u (0 , δ ˆ ψ ) − δψ (cid:17) H Γ == (cid:16) γ Γ u ( φ, ˆ ψ ) − E Γ ˆ ψ, γ Γ u (0 , δ ˆ ψ ) (cid:17) H Γ − (cid:16) γ Γ u ( φ, ˆ ψ ) − E Γ ˆ ψ, E Γ δ ˆ ψ (cid:17) H Γ ++ (cid:16) γ Γ E Σ ˆ u ( φ, ˆ ψ ) − E Γ ˆ ψ, γ Γ E Σ ˆ u (0 , δ ˆ ψ ) (cid:17) H Γ − (cid:16) γ Γ E Σ ˆ u ( φ, ˆ ψ ) − E Γ ˆ ψ, E Γ δ ˆ ψ (cid:17) H Γ == (cid:68) A ∗ p, A − Cδ ˆ ψ (cid:69) V (cid:48) ,V − (cid:68) E Γ ∗ Θ H Γ ( γ Γ u ( φ, ˆ ψ ) − E Γ ˆ ψ ) , δ ˆ ψ (cid:69) ˆ V (cid:48) , ˆ V ++ (cid:68) (cid:98) A ∗ ˆ p, (cid:98) A − (cid:98) Cδ ˆ ψ (cid:69) ˆ V (cid:48) , ˆ V − (cid:68) E Γ ∗ Θ H Γ ( E Γ ˆ u ( φ, ˆ ψ ) − E Γ ˆ ψ ) , δ ˆ ψ (cid:69) ˆ V (cid:48) , ˆ V == (cid:68) C ∗ p, δ ˆ ψ (cid:69) ˆ V (cid:48) , ˆ V − (cid:68) E Γ ∗ Θ H Γ ( γ Γ u ( δφ, ˆ ψ ) − E Γ ˆ ψ ) , δ ˆ ψ (cid:69) ˆ V (cid:48) , ˆ V ++ (cid:68) (cid:98) C ∗ ˆ p, δ ˆ ψ (cid:69) ˆ V (cid:48) , ˆ V − (cid:68) E Γ ∗ Θ H Γ ( E Γ ˆ u ( φ, ˆ ψ ) − E Γ ˆ ψ ) , δ ˆ ψ (cid:69) ˆ V (cid:48) , ˆ V == (cid:16) Θ − V ( C ∗ p + (cid:98) C ∗ ˆ p − E Γ ∗ Θ H Γ γ Γ u ( φ, ˆ ψ ) − E Γ ∗ Θ H Γ E Γ ˆ u ( φ, ˆ ψ ) + 2 E Γ ∗ Θ H Γ E Γ ˆ ψ ) , δ ˆ ψ (cid:17) ˆ V , which yield the thesis.Starting from the derivatives computed in Proposition 1 let us define the quan-tities δφ = Θ ˆ V ( B ∗ p − (cid:98) B ∗ ˆ p ) ∈ ˆ V (cid:48) (43) δ ˆ ψ = Θ − V ( C ∗ p + (cid:98) C ∗ ˆ p − E Γ ∗ Θ H Γ ( γ Γ u ( φ, ˆ ψ ) + E Γ ˆ u ( φ, ˆ ψ ) − E Γ ˆ ψ )) ∈ ˆ V . (44)Then the following proposition holds:
Proposition 2.
Given the variable χ = ( φ, ˆ ψ ) , let us increment it by a step ζδχ ,where δ X = ( δφ, δ ˆ ψ ) . The steepest descent method corresponds to the stepsize ζ = − (cid:0) δφ, δφ (cid:1) ˆ V (cid:48) + (cid:16) δ ˆ ψ, δ ˆ ψ (cid:17) ˆ V (cid:68) Bδφ + Cδ ˆ ψ, δp (cid:69) V (cid:48) ,V + (cid:68) − (cid:98) Bδφ + (cid:98) Cδ ˆ ψ, δ ˆ p (cid:69) ˆ V (cid:48) , ˆ V − (cid:68) E Γ ∗ Θ H Γ ( γ Γ δu + E Γ δ ˆ u − E Γ δ ˆ ψ ) , δ ˆ ψ (cid:69) ˆ V (cid:48) , ˆ V here δu = u ( δφ, δ ˆ ψ ) = A − ( Bδφ + Cδ ˆ ψ ) ∈ V,δ ˆ u = ˆ u ( δφ, δ ˆ ψ ) = (cid:98) A − ( − (cid:98) Bδφ + (cid:98) Cδ ˆ ψ ) ∈ ˆ V and δp ∈ V , δ ˆ p ∈ ˆ V are such that: A ∗ δp = γ ∗ Γ Θ H Γ ( γ Γ δu − E Γ δ ˆ ψ ) (cid:98) A ∗ δ ˆ p = E Γ ∗ Θ H Γ ( E Γ δ ˆ u − E Γ δ ˆ ψ ) Proof.
It is sufficient to set to zero the derivative ∂J ( χ + ζδχ ) ∂ζ . In order to get alighter notation let us set u = u ( φ, ˆ ψ ); δu = u ( δφ, δ ˆ ψ ); ˆ u = ˆ u ( φ, ˆ ψ ); δ ˆ u = ˆ u ( δφ, δ ˆ ψ ); J ( χ + ζδχ ) = J ( φ + ζδφ, ˆ ψ + ζδ ˆ ψ ) == 12 (cid:16) γ Γ u ( φ + ζδφ, ˆ ψ + ζδ ˆ ψ ) − ψ − ζδψ, γ Γ u ( φ + ζδφ, ˆ ψ + ζδ ˆ ψ ) − ψ − ζδψ (cid:17) H Γ + 12 (cid:16) γ Γ ˜ u ( φ + ζδφ, ˆ ψ + ζδ ˆ ψ ) − ψ − ζδψ, γ Γ ˜ u ( φ + ζδφ, ˆ ψ + ζδ ˆ ψ ) − ψ − ζδψ (cid:17) H Γ == 12 (cid:16) γ Γ u + ζγ Γ δu − E Γ ˆ ψ − ζ E Γ δ ˆ ψ, γ Γ u + ζγ Γ δu − E Γ ˆ ψ − ζ E Γ δ ˆ ψ (cid:17) H Γ ++ 12 (cid:16) γ Γ E Σ ˆ u + ζγ Γ E Σ δ ˆ u − E Γ ˆ ψ − ζ E Γ δ ˆ ψ, γ Γ E Γ ˆ u + ζγ Γ E Γ δ ˆ u − E Γ ˆ ψ − ζ E Γ δ ˆ ψ (cid:17) H Γ == J ( φ, ˆ ψ ) + ζ (cid:16) γ Γ u − E Γ ˆ ψ, γ Γ δu − E Γ δ ˆ ψ (cid:17) H Γ + ζ (cid:16) E Γ ˆ u − E Γ ˆ ψ, E Γ δ ˆ u − E Γ δ ˆ ψ (cid:17) H Γ ++ ζ (cid:16) γ Γ δu − E Γ δ ˆ ψ, γ Γ δu − E Γ δ ˆ ψ (cid:17) H Γ + ζ (cid:16) E Γ δ ˆ u − E Γ δ ˆ ψ, E Γ δ ˆ u − E Γ δ ˆ ψ (cid:17) H Γ ∂J ( χ + ζδχ ) ∂ζ = (cid:16) γ Γ u − E Γ ˆ ψ, γ Γ δu − E Γ δ ˆ ψ (cid:17) H Γ + (cid:16) E Γ ˆ u − E Γ ˆ ψ, E Γ δ ˆ u − E Γ δ ˆ ψ (cid:17) H Γ ++ ζ (cid:16) γ Γ δu − E Γ δ ˆ ψ, γ Γ δu − E Γ δ ˆ ψ (cid:17) H Γ + ζ (cid:16) E Γ δ ˆ u − E Γ δ ˆ ψ, E Γ δ ˆ u − E Γ δ ˆ ψ (cid:17) H Γ = 0 ⇒ ζ = − (cid:16) γ Γ u − E Γ ˆ ψ, γ Γ δu − E Γ δ ˆ ψ (cid:17) H Γ + (cid:16) E Γ ˆ u − E Γ ˆ ψ, E Γ δ ˆ u − E Γ δ ˆ ψ (cid:17) H Γ (cid:16) γ Γ δu − E Γ δ ˆ ψ, γ Γ δu − E Γ δ ˆ ψ (cid:17) H Γ + (cid:16) E Γ δ ˆ u − E Γ δ ˆ ψ, E Γ δ ˆ u − E Γ δ ˆ ψ (cid:17) H Γ ζ = − (cid:68) A ∗ p, A − ( Bδφ + Cδ ˆ ψ ) (cid:69) V (cid:48) ,V + (cid:68) (cid:98) A ∗ ˆ p, (cid:98) A − ( − (cid:98) Bδφ + (cid:98) Cδ ˆ ψ ) (cid:69) ˆ V (cid:48) , ˆ V − (cid:68) E Γ ∗ Θ H Γ ( γ Γ u + E Γ ˆ u − E Γ ˆ ψ ) , δ ˆ ψ (cid:69) ˆ V (cid:48) , ˆ V (cid:68) A − ( Bδφ + Cδ ˆ ψ ) , A ∗ δp (cid:69) V,V (cid:48) + (cid:68) (cid:98) A − ( − (cid:98) Bδφ + (cid:98) Cδ ˆ ψ ) , (cid:98) A ∗ δ ˆ p (cid:69) ˆ V , ˆ V (cid:48) − (cid:68) E Γ ∗ Θ H Γ ( γ Γ δu + E Γ δ ˆ u − E Γ δ ˆ ψ, δ ˆ ψ ) (cid:69) ˆ V (cid:48) , ˆ V == − (cid:68) B ∗ p, δφ (cid:69) ˆ V , ˆ V (cid:48) + (cid:68) C ∗ p, δ ˆ ψ (cid:69) ˆ V (cid:48) , ˆ V − (cid:68) (cid:98) B ∗ ˆ p, δφ (cid:69) ˆ V , ˆ V (cid:48) + (cid:68) (cid:98) C ∗ ˆ p, δ ˆ ψ (cid:69) ˆ V (cid:48) , ˆ V − (cid:68) E Γ ∗ Θ H Γ ( γ Γ u + E Γ ˆ u − E Γ ˆ ψ ) , δ ˆ ψ (cid:69) ˆ V (cid:48) , ˆ V (cid:68) Bδφ + Cδ ˆ ψ, δp (cid:69) V (cid:48) ,V + (cid:68) − (cid:98) Bδφ + (cid:98) Cδ ˆ ψ, δ ˆ p (cid:69) ˆ V (cid:48) , ˆ V − (cid:68) E Γ ∗ Θ H Γ ( γ Γ δu + E Γ δ ˆ u − E Γ δ ˆ ψ ) , δ ˆ ψ (cid:69) ˆ V (cid:48) , ˆ V == − (cid:16) δφ, δφ (cid:17) ˆ V (cid:48) + (cid:16) δ ˆ ψ, δ ˆ ψ (cid:17) ˆ V (cid:68) Bδφ + Cδ ˆ ψ, δp (cid:69) V (cid:48) ,V + (cid:68) − (cid:98) Bδφ + (cid:98) Cδ ˆ ψ, δ ˆ p (cid:69) ˆ V (cid:48) , ˆ V − (cid:68) E Γ ∗ Θ H Γ ( γ Γ δu + E Γ δ ˆ u − E Γ δ ˆ ψ ) , δ ˆ ψ (cid:69) ˆ V (cid:48) , ˆ V that yields the thesis.
5. Managing multiple cylinders and their intersections
The previous discussion can be readily adapted to the case of multiple, say I , small cylindrical inclusions Σ k , with lateral surface Γ k and centreline Λ i , i = 1 , . . . , I . Function spaces ˆ V , (cid:101) V , H Γ are introduced on each segment anddenoted by ˆ V i , (cid:101) V i , H Γ i , while trace operator γ Γ , extension operator E Γ and E Σ areeasily re-defined for each segment and denoted by γ Γ i , E Γ i and E Σ i , respectively.Operators (27)-(32) are re-written as: (cid:104) Au, v (cid:105) V (cid:48) ,V = ( K ∇ u, ∇ v ) V + α I (cid:88) i =1 ( | Γ i | ˇ u i , ˇ v i ) ˆ V i ∀ v ∈ V, ˇ v i ∈ ˆ V i : γ Γ i v = E Γ i ˇ v i , ∀ i = 1 , . . . , I (cid:68) (cid:98) A i ˆ u i , ˆ v i (cid:69) ˆ V (cid:48) i , ˆ V i = (cid:16) ˜ K i | Σ i | d ˆ u i ds , d ˆ v i ds (cid:17) ˆ V i + ˆ α ( | Γ i | ˆ u i , ˆ v i ) ˆ V i ∀ ˆ v i ∈ ˆ V i (cid:10) B i φ i , v i (cid:11) V (cid:48) i ,V i = (cid:10) | Γ i | φ i , ˇ v i (cid:11) ˆ V (cid:48) i , ˆ V i ∀ v ∈ V : γ Γ i v = E Γ i ˇ v i , ˇ v i ∈ ˆ V i (cid:68) (cid:98) B i φ i , ˆ v (cid:69) ˆ V (cid:48) i , ˆ V i = (cid:10) | Γ i | φ i , ˆ v (cid:11) ˆ V (cid:48) i , ˆ V i ∀ ˆ v ∈ ˆ V i (cid:68) C i ˆ ψ i , v (cid:69) V (cid:48) ,V = α ( | Γ i | ˆ ψ, ˇ v ) ˆ V i ∀ v ∈ V, ˇ v i ∈ ˆ V i : γ Γ i v = E Γ i ˇ v (cid:68) (cid:98) C i ˆ ψ i , ˆ v (cid:69) ˆ V (cid:48) i , ˆ V i = ˆ α ( | Γ i | ˆ ψ i , ˆ v ) ˆ V i ˆ v ∈ ˆ V i . Problem (35)-(36) can be finally re-written for an arbitrary set of I centrelinesas: Au − B i φ i − C i ˆ ψ i = F (45) (cid:98) A i ˆ u i + (cid:98) B i φ i − (cid:98) C i ˆ ψ i = G i , (46)11ith the definition of G i following from the definition of an operator in theform of (34) for each segment. in which we consider different 1D variables ˆ u i , φ i , ˆ ψ i on the different 1D domains. In case cylindrical inclusions withintersecting centrelines, i.e. such that ¯Λ i ∪ ¯Λ j (cid:54) = ∅ , i, j ∈ [1 , . . . , I ] , we canstill adopt formulation (45)-(46) by splitting the intersecting centrelines intonon intersecting sub-segments and then adding continuity conditions at theintersection points.Finally, the cost functional is re-written as: J = I (cid:88) i =1 J i := I (cid:88) k =1 (cid:16) || γ Γ i u ( φ i , ˆ ψ i ) − ψ i || H Γ + || γ Γ ˜ u i ( φ i , ˆ ψ i ) − ψ i || H Γ (cid:17) (47)where ˜ u i = E Σ i ˆ u i and ψ i = E Γ i ˆ ψ i .
6. Discrete matrix formulation
Here the discrete matrix form of problem (37)-(38) is presented. The 3D-1D coupling is trivial in the discrete approximation spaces, given the regularityproperties of the function spaces commonly used for discretization. Nonetheless,the present approach has the advantage of having a well posed mathematicalformulation, and it even allows the use of non conforming meshes at the inter-faces of the subdomains. Indeed, thanks to the optimization framework, it ispossible to use completely independent meshes for the various domains and alsofor the interface variables, without any theoretical constraint on mesh sizes.For the sake of generality we will consider, from the beginning, the presenceof multiple segments crossing domain Ω . Let us consider I segments of differentlength and orientation, defined ad Λ i = { λ i ( s ) , s ∈ (0 , S i ) } , i = 1 , ..., I . Let usconsider a tetrahedral mesh T of domain Ω , and let us define, on this mesh,Lagrangian finite element basis functions { ϕ k } Nk =1 , such that U = (cid:80) Nk =1 U k ϕ k is the discrete approximation of pressure u . Let us then build three partitions ofeach segment Λ i , named ˆ T i , τ φi and τ ψi , defined independently from each otherand from T . Let us further define the basis functions { ˆ ϕ i,k } ˆ N i k =1 on ˆ T i , { θ i,k } N φi k =1 on τ φi and { η i,k } N ψi k =1 on τ ψi , with ˆ N i , N φi and N ψi denoting the number of DOFsof the discrete approximations of the variables ˆ u i , φ i and ˆ ψ i respectively, havingset: ˆ U i = ˆ N i (cid:88) k =1 ˆ U i,k ˆ ϕ i,k , Φ i = N φi (cid:88) k =1 Φ i,k θ i,k , Ψ i = N ψi (cid:88) k =1 Ψ i,k η i,k We then define the following matrices: A ∈ R N × N s.t. ( A ) kl = (cid:90) Ω K ∇ ϕ k ∇ ϕ l dω + α I (cid:88) i =1 (cid:90) Λ i | Γ( s i ) | ϕ k | Λ i ϕ l | Λ i ds ˆ A i ∈ R ˆ N i × ˆ N i s.t. ( ˆ A i ) kl = (cid:90) Λ i ˜ K i | Σ( s i ) | d ˆ ϕ i,k ds d ˆ ϕ i,l ds ds + ˆ α (cid:90) Λ i | Γ( s i ) | ˆ ϕ i,k ˆ ϕ i,l ds i ∈ R N × N φi s.t. ( B i ) kl = (cid:90) Λ i | Γ( s i ) | ϕ k | Λ i θ i,l ds ˆ B i ∈ R ˆ N i × N φi s.t. ( ˆ B i ) kl = (cid:90) Λ i | Γ( s i ) | ˆ ϕ i,k θ i,l ds C α i ∈ R N × N ψi s.t. ( C αi ) kl = α (cid:90) Λ i | Γ( s i ) | ϕ k | Λ i η i,l ds ˆ C α i ∈ R ˆ N i × N ψi s.t. ( ˆ C iα ) kl = ˆ α (cid:90) Λ i | Γ( s i ) | ˆ ϕ i,k η i,l ds, and the vectors f ∈ R N s.t. f k = (cid:90) Ω f ϕ k dωg i ∈ R ˆ N i s.t. ( g i ) k = (cid:90) Λ i | Σ( s i ) | g ˆ ϕ i,k ds. Setting ˆ N = (cid:80) I i =1 ˆ N i , N ψ = (cid:80) I i =1 N ψi and N φ = (cid:80) I i =1 N φi , we can group thematrices as follows for all the segments in the domain: B = [ B , B , ..., B I ] ∈ R N × N φ ˆ B = diag (cid:16) ˆ B , ..., ˆ B I (cid:17) ∈ R ˆ N × N φ C α = [ C α , C α , ..., C α I ] ∈ R N × N ψ ˆ C α = diag (cid:16) ˆ C α , ..., ˆ C α I (cid:17) ∈ R ˆ N × N ψ and, for non intersecting segments, we have: ˆ A = diag (cid:16) ˆ A , ..., ˆ A I (cid:17) ∈ R ˆ N × ˆ N , whereas, for groups of intersecting segments, we proceed as described in Sec-tion 5 and we enforce continuity through Lagrange multipliers. For each con-nected group of segments we thus have: ˆ A (cid:63)ζ = (cid:34) diag (cid:16) ˆ A ζ , ..., ˆ A ζ n (cid:17) Q T Q (cid:35) where matrix Q simply equates the DOFs at the extrema of connected sub-segments. Matrices ˆ A (cid:63)ζ , for ζ spanning the whole number of connected groupsof segments, are assembled block diagonally to form matrix ˆ A . Please note that,for disconnected segments, each matrix ˆ A (cid:63)ζ coincides with matrix ˆ A ζ . Finallywe can write A U − B Φ − C α Ψ = f (48) ˆ A ˆ U + ˆ B Φ − ˆ C α Ψ = g (49)with ˆ U = (cid:104) ˆ U T , ..., ˆ U T I (cid:105) T ∈ R ˆ N ; g = [ g T , g T , ..., g T I ] T ∈ R ˆ N Φ = (cid:2) Φ T , ..., Φ T I (cid:3) T ∈ R N φ ; Ψ = (cid:2) Ψ T , ..., Ψ T I (cid:3) T ∈ R N ψ .
13n order to get a more compact form of the previous equations, let us set W = ( U, ˆ U ) and A = (cid:20) A ˆ A (cid:21) , B = (cid:20) B − ˆ B (cid:21) , C α = (cid:20) C α ˆ C α (cid:21) F = (cid:20) fg (cid:21) , (50)so that the discrete constraint equations become: A W − B Φ + C α Ψ = F . (51)Concerning the cost functional in (26), first we define matrices G i ∈ R N × N s.t. ( G i ) kl = (cid:90) Λ i ϕ k | Λ i ϕ l | Λ i ds ˆ G i ∈ R ˆ N i × ˆ N i s.t. ( ˆ G i ) kl = (cid:90) Λ i ˆ ϕ i,k ˆ ϕ i,l ds G ψi ∈ R N ψi × N ψi s.t. ( G ψi ) kl = (cid:90) Λ i η i,k η i,l ds C i ∈ R N × N ψi s.t. ( C i ) kl = (cid:90) Λ i ϕ k | Λ i η i,l ds ˆ C i ∈ R ˆ N i × N ψi s.t. ( ˆ C i ) kl = (cid:90) Λ i ˆ ϕ i,k η i,l ds and then G = I (cid:88) i =1 G i ∈ R N × N ˆ G = diag (cid:16) ˆ G T , ..., ˆ G T I (cid:17) ∈ R ˆ N × ˆ N G = (cid:20) G ˆ G (cid:21) (52) G ψ = diag (cid:16) G ψ , ..., G ψ I (cid:17) ∈ R N ψ × N ψ C = [ C , C , ..., C I ] ∈ R N × N ψ ˆ C = diag (cid:16) ˆ C , ..., ˆ C I (cid:17) ∈ R ˆ N × N ψ C = (cid:20) C ˆ C (cid:21) . The discrete cost functional then reads: ˜ J = 12 (cid:16) U T G U − U T C Ψ − Ψ T C T U + ˆ U T ˆ G ˆ U − ˆ U T ˆ C Ψ − Ψ T ˆ C T ˆ U + 2Ψ T G ψ Ψ (cid:17) == 12 (cid:16) W T G W − W T C Ψ − Ψ T C T W + 2Ψ T G ψ Ψ (cid:17) . (53)The discrete matrix formulation of the 3D-1D problem finally takes the form: min (Φ , Ψ) ˜ J (Φ , Ψ) subject to (51) . (54)First order optimality conditions for the above problem correspond to the saddle-point system: 14 = G − C A T B T − C T G ψ ( − C α ) T A B − C α (55) S W ΦΨ − P = F (56) Proposition 3.
Matrix S in (55) is non-singular and the unique solution of (56) is equivalent to the solution of the optimization problem (54) . The proof of Proposition 3 derives from classical arguments of quadraticprogramming once the following lemma is proven:
Lemma 1.
Let matrix A (cid:63) be as A (cid:63) = (cid:2) A B − C α (cid:3) and let G (cid:63) be defined as G (cid:63) = G − C − C T G ψ Then matrix A (cid:63) is full rank and matrix G (cid:63) is symmetric positive definite on ker( A (cid:63) ) .Proof. The proof is adapted from the one provided in [18], we report here the keysteps. Matrix A is full rank and matrix G (cid:63) is symmetric positive semi-definiteby construction. We thus need to show that ker ( G (cid:63) ) ∩ ker ( A (cid:63) ) = { } . Let usconsider the canonical basis for R N Φ + N Ψ and let e k denote the k -th element ofsuch basis, k = 1 , . . . , N Φ + N Ψ . Let z k ∈ ker ( A (cid:63) ) be defined as: z k = (cid:20) A − (cid:2) B − C α (cid:3) e k e k (cid:21) . Let us assume that ≤ k ≤ N Φ , thus corresponding to a non null value of thevariable Φ on one segment. This in turn gives a non null value U and ˆ U on thetraces and thus a non null value of the functional, or z Tk G (cid:63) z k > . If instead N Φ + 1 ≤ k ≤ N Ψ , this corresponds to a non-null variable Ψ := e k . Since thesolution to (56) is the same for every value of α and ˆ α , included α = ˆ α = 0 (theconsistent terms depending on α and ˆ α are only required for the independentresolution on the sub-domains), we choose here α = ˆ α = 0 , so that: z k = (cid:20) A − (cid:2) B (cid:3) e k e k (cid:21) := (cid:20) e k (cid:21) U , ˆ U are null, and therefore we can conclude that z Tk G z k > also inthis case (see [18] for the proof with α, ˆ α > ). Summarizing we have shownthat z k (cid:54)∈ ker ( G (cid:63) ) for any k = 1 , . . . , N Φ + N Ψ . The vector space ker ( A (cid:63) ) = span { z , . . . , z N Φ + N Ψ } is a subspace of Im ( G (cid:63) ) , and ker ( G (cid:63) ) ∩ ker ( A (cid:63) ) = { } .System (56) can be used to obtain a numerical solution to the problem. Forvery large problems instead, the above system is likely to become ill-conditioned,and thus an alternative resolution strategy is proposed, as described below. Byformally replacing W = A − ( B Φ − C α Ψ + F ) in the functional (53), we obtain J (cid:63) (Φ , Ψ) = 12 (cid:16) ( A − B Φ + A − C α Ψ + A − F ) T G ( A − B Φ + A − C α Ψ + A − F )+ − ( A − B Φ + A − C α Ψ + A − F ) T C Ψ+ − Ψ T C T ( A − B Φ + A − C α Ψ + A − F ) (cid:17) == 12 [Φ T Ψ T ] B T A − T GA − B B T A − T GA − C α + − B T A − T C ( C α ) T A − T GA − B + − C T A − B ( C α ) T A − T GA − C α + − C T A − T C α + − ( C α ) T A − C +2 G ψ ΦΨ ++ F T (cid:2) A − T GA − B A − T GA − C α − A − T C (cid:3) (cid:20) ΦΨ (cid:21) ++ 12 (cid:16) F T A − T GA − F (cid:17) == 12 (cid:0) X T M X + 2 d T X + q (cid:1) . (57)Matrix M is symmetric positive definite as follows from the equivalence of thisformulation with the previous saddle point system (56), and thus the minimiza-tion of the unconstrained problem (57) can be performed via a gradient basedscheme. It is however to remark that the computation of gradient direction atpoint X (cid:93) , i.e. ∇ J (cid:63) ( X (cid:93) ) = M X (cid:93) + d , can be performed in a matrix free formatand involves the independent factorization of the 1D matrices ˆ A i , i = 1 , . . . , L and of the 3D elliptic matrix A , which are all non singular as long as α, ˆ α > .The analysis of this solving strategy and of its potential for parallel computingis deferred to a forthcoming work.
7. Numerical results
In this section we propose some numerical test to validate the proposed ap-proach and to show its applicability to the problem of interest. Three numericaltests are proposed. A first problem called
Test Problem 1 (TP1) takes intoaccount a single cylindrical inclusion and has a smooth analytical solution, thus16 nduced ˆ δ u =1 ˆ δ u =0 . Figure 1: Highlight on the tetrahedra intersected by a segment and consequently inducedmesh; on the right equispaced partitions of the segment for ˆ δ u = 1 and for ˆ δ u = 0 . . allowing to evaluate convergence trends for the error. The second test, called Test Problem 2 (TP2), takes into account a different problem with no knownanalytical solution on a similar geometry. In this case, the obtained solution iscompared to a 3D-3D simulation with standard conforming finite elements, anddifferent values of the coefficient ˜ K are considered. Finally an example withmultiple intersecting inclusions is proposed, to test the behavior of the methodin more general settings. In this case, a qualitative evaluation on the behaviorof the numerical solution is proposed, along with a quantitative evaluation of aproposed error indicator.All the simulations are performed using finite elements on 3D and 1D non-conforming meshes, independently generated on the sub-domains. A mesh pa-rameter h is used to denote the maximum diameter of the tetrahedra for the 3Dmesh of Ω , whereas the refinement level of the 1D meshes ˆ T i , τ φi , τ ψi , i = 1 , ... I ,is provided in relation to the mesh-size of the 1D mesh induced on the segments Λ i by the tetrahedral mesh, i.e. the 1D mesh given by the intersections of Λ i with the tetrahedra in T , see Figure 1. This is done in order to better high-light the relative sizes of the various meshes. In particular the adimensionalnumber ˆ δ u,i denotes the ratio between the number of elements of the mesh in ˆ T i with respect to the number of elements of the induced mesh on Λ i , whereas δ φ,i and δ ψ,i the ratio between the number of elements in τ φi , τ ψi , respectively,and the induced mesh. Figure 1 shows the induced mesh, in the middle, andtwo 1D meshes, e.g. for ˆ T corresponding to values of ˆ δ u = 1 and ˆ δ u = 0 . and equally-spaced nodes. In the simulations, for simplicity, we will always useequally-spaced nodes for the 1D meshes and unique different values of ˆ δ u , δ φ and δ ψ for the various segments, thus dropping, in the following, the referenceto segment index for these parameters. Linear Lagrangian finite elements are17 Σ Figure 2: TP1: on the left, section view of the starting 3D-3D domain for the TestProblem1experiment; on the right solution obtained inside the cube and on the segment for h = 0 . , ˆ δ u = 1 , δ φ = 0 . and δ ψ = 0 . . used on T , ˆ T i , τ ψi , whereas piece-wise constant basis functions are used to de-scribe variables φ i on τ φi , i = 1 , ... I . All numerical tests are performed using α = ˆ α = 1 , even if it is to remark that the value of such parameters has noimpact on the solution, and, as long as formulation (56) is used, α = ˆ α = 0 could have been also chosen. Let us consider a cube Ω of edge l inscribed in a cylinder of radius ˆ R = l √ centered in the axis origin and a cylinder Σ of radius ˇ R < ˆ R and height h = l whose centreline Λ lies on the z axis (see Figure 2). Let us denote by ∂ Ω l , ∂ Ω + and ∂ Ω − respectively the lateral, the top and the bottom faces of the cube.Given a, b, c, k , k ∈ R , let us consider a problem in the form of (20)-(23),obtained by reducing Σ to its centerline, with K = ˜ K = 1 and f = − b (cid:112) x + y − a, g = 0 . The problem is completed with the appropriate boundary condition to havethe exact solution given by: u ex ( x, y, z ) = a ( x + y ) + b (cid:112) x + y + c in Ω (58) ˆ u ex ( x, y, z ) = k on Λ (59)with a = k − k ( ˆ R − ˇ R ) , b = − R ( k − k )( ˆ R − ˇ R ) , c = k + ( k − k ) ˇ R ( ˆ R − ˇ R ) This reduced problem corresponds to an equi-dimensional problem satisfyingour modeling assumptions and having a constant solution equal to k insidethe cylinder. Further, the flux through the interface is zero, as the solution is C in the whole domain. Results are obtained considering a cube of edge l = 2 igure 3: TP1: trend of the L and H -norms of the relative errors under mesh refinement.On the left: error committed on the cube with respect to (58); on the right: error committedon the segment with respect to (59). Other parameters: ˆ δ u = 1 , δ ψ = δ φ = 0 . . ( ˆ R = √ ) and choosing ˇ R = 0 . , k = 0 . and k = 5 . Homogeneous Neumannboundary conditions are imposed on ∂ Ω + and ∂ Ω − , whereas Dirichlet boundaryconditions are imposed on ∂ Ω l . Dirichlet boundary conditions equal to k areimposed on segment endpoints. Figure 2 on the right shows the approximatedsolutions U , ˆ U obtained inside the cube and on the segment for h = 0 . and ˆ δ u = 1 , corresponding to N = 3715 DOFs in the cube and ˆ N = 57 DOFs on thesegment. The other parameters are δ φ = 0 . and δ ψ = 0 . . Convergence curvesof the error can be computed and, given the regularity of the solution, optimalconvergence trends are expected for the used finite element approximation. Letus introduce the errors E L , E H , for the 3D problem and (cid:98) E L and (cid:98) E H for the1D problem, defined as follows: E L = || u ex − U || L (Ω) || u ex || L (Ω) , E H = || u ex − U || H (Ω) || u ex || H (Ω) , (cid:98) E L = || ˆ u ex − ˆ U || L (Λ) || ˆ u ex || L (Λ) , (cid:98) E H = || ˆ u ex − ˆ U || H (Λ) || ˆ u ex || H (Λ) . Figure 3 displays the convergence trends for the above quantities against meshrefinement. Four meshes are considered, characterized by mesh parameters h =0 . , . , . , . respectively, corresponding to N = 229 , , , DOFs and ˆ N = 20 , , , DOFs ( ˆ δ u = 1 , δ φ = δ ψ = 0 . ), confirming the ex-pected behaviors.For this simple problem it is possible to compute the condition number ofthe KKT optimality conditions (55), and analyze how it is affected by differentchoices of the meshsize of the 1D meshes. Figure 4, on the left, reports thecondition number of the KKT system matrix as δ φ varies between . and ,for five different values of ˆ δ u between . and . , being instead δ ψ = 0 . fixed.We can see that δ φ has a relatively small impact on the conditioning of thesystem as long as δ φ < ˆ δ u is chosen, otherwise a large rapid increase is observedas δ φ ≥ ˆ δ u grows. It is therefore advisable to choose a quite coarse mesh forvariable Φ with respect to the mesh for ˆ U , even if no theoretical constraints19 igure 4: TP1: trend of the conditioning of the KKT system under the variation of the 1Dmesh parameters. On the left variable δ φ and different values of ˆ δ u , while δ ψ = 0 . . On theright variable δ ψ and δ φ = 0 . . In both cases h = 0 . . emerged in the analysis. Figure 4, on the right, shows again the conditioning ofthe KKT system matrix as δ ψ varies between . and , for the same five valuesof ˆ δ u between . and . , keeping this time δ φ = 0 . fixed. It can be seen that,in this case, the conditioning is almost independent of δ ψ , for all the consideredvalues of ˆ δ u . It is however to remark that saddle point matrices as the one in(55) are typically ill conditioned. The use of a resolution strategy based on agradient based scheme for the minimization of the unconstrained functional isexpected to result in a problem with a mitigated condition number, actuallycoinciding with the application of a null-space based preconditioning technique[22]. The second example is set on a domain equal to the one of TP1. We considerthree different problems defined as in (20)-(23), characterized by three differentvalues of the coefficient ˜ K , equal to , and respectively, whereas K = 1 , f = 1 , g = 0 are fixed for all the problems. Even the boundary conditions areshared: being ∂ Ω + , ∂ Ω − , and ∂ Ω l defined as previously, homogeneous Dirichletboundary conditions are prescribed on ∂ Ω + , ∂ Ω − and at segment endpoints,while homogeneous Neumann boundary conditions are set on ∂ Ω l .The accuracy of the solution is evaluated by means of a comparison withan equi-dimensional problem having a cylindrical inclusion of radius . , withcentreline coinciding with the 1D domain of the reduced problem. Dirichlethomogeneous boundary conditions are prescribed on the top and bottom facesof the 3D domains and homogeneous Neumann boundary conditions are set onthe outer surface. A unitary forcing term is prescribed in the 3D domain outsideof the cylindrical inclusion, where, instead a null forcing is set. As ˜ K grows, wemove from a problem with a smooth solution to a problem with a jump of thegradient across the interfaces between the bulk 3D domain and the inclusion.The equi-dimensional problem is solved on a fine mesh, refined around theinclusion in order to match the geometry. The lateral surface of the cylindrical20 igure 5: TP2: Top: mesh used for the reference equi-dimensional problem, conforming to thecylindrical inclusion. Bottom: adapted non conforming mesh for the 3D-1D reduced problemwith the proposed approachFigure 6: TP2: Uniformly refined mesh with h = 0 . . igure 7: TP2: comparison of the results obtained along the centerline of the cylinder in the3D-3D conforming setting and by using the 3D-1D reduced model. Other parameters for the3D-1D setting: ˆ δ u = 1 and δ φ = δ ψ = 0 . .Figure 8: TP2: comparison of the solution obtained on the adapted mesh for the 3D-1Dproblem with the reference solution on a plane parallel to y − z and containing the centrelineof the inclusion. ˜ K = 10 . . million elements and DOFs.The reduced 3D-1D problem is solved on four different meshes: first a nonconforming mesh, slightly refined close to the inclusion area, is considered,termed adapted mesh and having . × elements and DOFs. Thismesh is thus much coarser than the reference mesh. It is shown in Figure 5, atthe bottom, along with a zooming of the zone around the 1D domain to high-light the non conformity. Further, three uniformly refined meshes are consid-ered, with mesh parameters h = 0 . , . , . , respectively, correspondingto N = 1287 , , DOFs. The intermediate uniform mesh is shown inFigure 6.The solution on the centreline of the inclusion obtained on the various consid-ered meshes are compared to the reference solution on the centreline in Figure 7,for ˜ K = 1 on the top left, ˜ K = 100 on the top right, and for ˜ K = 10 on thebottom. We can clearly see that, as long as the jump in the coefficient betweenthe bulk domain and the inclusion is relatively small, the proposed approachcorrectly reproduces the solution on all the considered meshes. Instead, forlarge jumps, as it is for ˜ K = 10 , the solution on the uniformly refined meshesare less accurate, whereas, the use of a slightly adapted mesh, even if still nonconforming, is capable of producing a solution in very good agreement with thereference. The proposed approach can be thus of help in mitigating the over-head in mesh generation and to reduce problem size. A comparison betweenthe reference solution and the solution of the reduced 3D-1D problem on theadapted mesh and with ˜ K = 10 is finally shown in Figure 8, on a slice parallelto the y − z plane and containing the centreline. The plot of the two solutionsmatch well. As for the previous cases, let us consider a cube of edge l = 2 centered in theaxes origin. Let us then consider a set of inclusions of radius ˇ R = 10 − , whosecenterlines intersect in points. We impose homogeneous Dirichlet boundaryconditions on all the faces of the cube and at the dead ends of the networkintersecting cube top and bottom faces, as shown in Figure 9. HomogeneousNeumann boundary conditions are imposed at segment endpoints lying insidethe cube. We consider a problem in the same form of (45)-(46), with i =1 , ..., I = 19 , spanning the segments, that form a unique connected component,as discussed in Section 6. Further, we consider K = 1 , f = 0 and ˜ K i = 100 , g i = 3 . e − ∀ i = 1 , ..., .We refer again to Figure 9 for the solution ˆ U obtained in the segment networkfor h = 0 . and ˆ δ u = 1 , corresponding to N = 12873 and ˆ N = 309 DOFs.The other parameters are δ φ = δ ψ = 0 . . Figure 10, on the left, shows instead,for the same choice of parameters, the solution U obtained inside the cube, onthree different slices, all parallel to the x − y plane and located at z = − . ,23 igure 9: MI: Solution obtained on the centerlines of the inclusions for h = 0 . , ˆ δ u = 1 and δ φ = δ ψ = 0 . . Homogeneous Dirichlet boundary conditions are imposed at the pointsmarked in blue.Figure 10: MI: On the left, solution obtained inside the cube on three different slices parallelto the x − y plane and located at z = − . , z = 0 and z = 0 . . On the right focus on thebehavior of the − K ∇ U vector field. The chosen mesh parameters are h = 0 . , ˆ δ u = 1 , δ φ = δ ψ = 0 . . z = 0 and z = 0 . . On the right of Figure 10 the behavior of the vector field − K ∇ U is shown.A quantitative evaluation of the numerical solution is provided through anerror indicator, denoted by ∆ L u , measuring how well the continuity condition,enforced through the minimization of the functional (53), is satisfied. We thusintroduce the quantity: ∆ L u = (cid:113)(cid:80) I i =1 || U | Λ i − ˆ U i || L (Λ i ) | max( U, ˆ U ) |√ l tot (60)resulting in a non dimensional number, with l tot denoting the total length ofthe segments in the domain.The trend of ∆ L u against δ φ and δ ψ , both ranging between . and , is24 igure 11: MI: Value of ∆ L u (60) under the variation of the mesh parameters. On theleft variable δ φ and δ ψ , h = 0 . ; on the right, variable ˆ δ u and four different values of h , δ φ = δ ψ = 0 . . shown in Figure 11 on the left, still considering h = 0 . and ˆ δ u = 1 . Asexpected, the error indicator decreases as the two parameters increase. Theimpact of δ φ appears to be stronger: for values of δ ψ close to , almost twoorders of magnitude are swept by ∆ L u as δ φ varies. The impact of δ ψ on thecontinuity condition appears to be weaker, even if it can be seen that it becomesmore relevant for high values of δ φ , with almost one order of magnitude sweptby the error indicator as δ ψ increases. Figure 11, on the right, shows insteadthe trend of ∆ L u against ˆ δ u , ranging between 0.6 and 2, for four values ofthe mesh size h , namely h = 0 . , . , . , . , corresponding to N =126 , , , . The other parameters are δ φ = δ ψ = 0 . . We can see thatthe continuity error indicator is only marginally affected by the value of ˆ δ u ,whereas, it can be arbitrarily reduced by mesh refinement. It should be noted,however, that, for a fixed value of ˆ δ u , a refinement of the 3D mesh also impliesa refinement of the 1D mesh for ˆ u , whereas, changes in ˆ δ u , at fixed h only refinethe mesh of ˆ u , leaving the 3D mesh unchanged.
8. Conclusions
A novel approach for 3D-1D coupled problems has been proposed. Themethod derives from a mathematical formulation in proper functional spacesthat allows the definition of a well posed trace operator from functions in thethree dimensional space to one dimensional manifolds. The 1D problems aredecoupled from the problem on the bulk 3D domain and two interface vari-ables are introduced in this domain decomposition process, thus resulting in athree field formulation of the original problem. A cost functional is introducedand minimized to impose matching conditions at the interfaces. The methodallows to enforce continuity conditions and flux balance at the interfaces be-tween sub-problems on non-conforming meshes, thus strongly alleviating themesh generation process. Indeed meshes on the various sub-domains can beindependently generated. Numerical results on two simple test problem and on25 more complex configuration show the viability of the proposed approach andare used to analyze the effect of method parameters on the condition numberof the discrete problem and on solution accuracy.A formulation suitable for efficient resolution through iterative gradient-based schemes is also envisaged and should be further investigated to allowsimulation on problems of high geometrical complexity.
ReferencesReferences [1] D. Notaro, L. Cattaneo, L. Formaggia, A. Scotti, P. Zunino, A Mixed FiniteElement Method for Modeling the Fluid Exchange Between Microcircula-tion and Tissue Interstitium, Springer International Publishing, 2016, pp.3–25. doi:10.1007/978-3-319-41246-7\_1 .[2] T. Köppl, E. Vidotto, B. Wohlmuth, A 3d-1d coupled blood flow andoxygen transport model to generate microvascular networks, InternationalJournal for Numerical Methods in Biomedical Engineering 36 (10) (2020)e3386. doi:10.1002/cnm.3386 .[3] N. Schröder, M. Javaux, J. Vanderborght, B. Steffen, H. Vereecken, Effectof root water and solute uptake on apparent soil dispersivity: A simulationstudy, Vadose Zone Journal 11 (3) (2012) vzj2012.0009. doi:https://doi.org/10.2136/vzj2012.0009 .[4] T. Koch, K. Heck, N. Schröder, H. Class, R. Helmig, A new simulationframework for soil–root interaction, evaporation, root growth, and so-lute transport, Vadose Zone Journal 17 (1) (2018) 170210. doi:10.2136/vzj2017.12.0210 .[5] I. G. Gjerde, K. Kumar, J. M. Nordbotten, Well modelling by means ofcoupled 1d-3d flow models, in: ECMOR XVI - 16th European Conferenceon the Mathematics of Oil Recovery, 2018.[6] I. Gjerde, K. Kumar, J. Nordbotten, A singularity removal method forcoupled 1d–3d flow models, Comput Geosci 24 (2020) 443–457. doi:10.1007/s10596-019-09899-4 .[7] D. Cerroni, F. Laurino, P. Zunino, Mathematical analysis, finite elementapproximation and numerical solvers for the interaction of 3d reservoirswith 1d wells, GEM - International Journal on Geomathematics 10 (1)(2019).[8] I. Steinbrecher, M. Mayr, M. Grill, J. Kremheller, C. Meier, A. Popp,A mortar-type finite element approach for embedding 1d beams into3d solid volumes, Comput Mech 66 (2020) 1377–1398. doi:10.1007/s00466-020-01907-0 . 269] A. Llau, L. Jason, F. Dufour, J. Baroth, Finite element modelling of 1d steelcomponents in reinforced and prestressed concrete structures, EngineeringStructures 127 (2016) 769–783. doi:10.1016/j.engstruct.2016.09.023 .[10] C. D’Angelo, Finite element approximation of elliptic problems withdirac measure terms in weighted spaces: applications to one- and three-dimensional coupled problems, SIAM J. Numer. Anal. 50 (1) (2012) 194 –215.[11] C. D’Angelo, A. Quarteroni, On the coupling of 1d and 3d diffusion-reactionequations. application to tissue perfusion problems, Math. Models MethodsAppl. Sci. 18 (2008) 1481 – 1504.[12] A. Ern, J. Guermond, Theory and Practice of Finite Elements, Vol. 159,Appl. Mat. Sci, Springer-Verlag, New York, 2004.[13] A.-K. Tornberg, B. Engquist, Numerical approximations of singular sourceterms in differential equations, Journal of Computational Physics 200 (2)(2004) 462–488. doi:10.1016/j.jcp.2004.04.011 .[14] Gjerde, Ingeborg G., Kumar, Kundan, Nordbotten, Jan M., Wohlmuth,Barbara, Splitting method for elliptic equations with line sources, ESAIM:M2AN 53 (5) (2019) 1715–1739. doi:10.1051/m2an/2019027 .[15] T. Köppl, E. Vidotto, B. Wohlmuth, P. Zunino, Mathematical modeling,analysis and numerical approximation of second-order elliptic problemswith inclusions, Mathematical Models and Methods in Applied Sciences28 (05) (2018) 953–978. doi:10.1142/S0218202518500252 .[16] F. Laurino, P. Zunino, Derivation and analysis of coupled pdes on manifoldswith high dimensionality gap arising from topological model reduction.,ESAIM: M2AN 53 (6) (2019) 2047 – 2080.[17] F. Brezzi, L. Marini, A three-field domain decomposition method, Contem-porary Mathematics 157 (1994).[18] S. Berrone, D. Grappein, S. Pieraccini, S.Scialò, A three-field based opti-mization formulation for flow simulations in networks of fractures on nonconforming meshes, accepted for publication on SIAM J. Sci. Comput.(2019). arXiv:1912.09744 .[19] S. Berrone, S. Pieraccini, S. Scialò, An optimization approach for largescale simulations of discrete fracture network flows, J. Comput. Phys. 256(2014) 838–853. doi:10.1016/j.jcp.2013.09.028 .[20] S. Berrone, S. Scialò, F. Vicini, Parallel meshing, discretization and compu-tation of flow in massive Discrete Fracture Networks, SIAM J. Sci. Comput.41 (4) (2019) C317–C338. doi:10.1137/18M1228736 .2721] S. Berrone, A. D’Auria, S. Scialò, An optimization approach for flow sim-ulations in poro-fractured media with complex geometries, Comput Geosci(2021). doi:10.1007/s10596-020-10029-8 .[22] J. Pestana, T. Rees, Null-space preconditioners for saddle point systems,SIAM Journal on Matrix Analysis and Applications 37 (3) (2016) 1103–1128. doi:10.1137/15M1021349doi:10.1137/15M1021349