# Linking discrete and continuum diffusion models: Well-posedness and stable
finite element discretizations

Christina Schenk1, David Portillo2,1, Ignacio Romero2,1 Corresponding author:
<EMAIL_ADDRESS>
(1IMDEA Materials Institute, Eric Kandel 2, Tecnogetafe, 28906 Madrid, Spain
2Universidad Politécnica de Madrid, José Gutiérrez Abascal, 2, 28006 Madrid,
Spain)

###### Abstract

In the context of mathematical modeling, it is sometimes convenient to
integrate models of different nature. These types of combinations, however,
might entail difficulties even when individual models are well-understood,
particularly in relation to the well-posedness of the ensemble. In this
article, we focus on combining two classes of dissimilar diffusive models: the
first one defined over a continuum and the second one based on discrete
equations that connect average values of the solution over disjoint
subdomains. For stationary problems, we show unconditional stability of the
linked problems and then the stability and convergence of its discretized
counterpart when mixed finite elements are used to approximate the model on
the continuum. The theoretical results are highlighted with numerical examples
illustrating the effects of linking diffusive models. As a side result, we
show that the methods introduced in this article can be used to infer the
solution of diffusive problems with incomplete data.

Keywords: Mixed finite elements; inf-sup condition; multi-scale; network
model; diffusion problems; stability

##  1 Introduction

Diffusion problems of interest in Applied Math and Engineering can be studied
with discrete models as well as partial differential equations (PDEs). The
first approach is naturally simpler than the second one since it sidesteps the
difficulties that result from the spatial description of the solution fields.
In addition, the solution of these simple models can be approximated very
efficiently at the expense of all spatial details, in contrast with the
approximation of PDEs. The latter, in fact, invariably requires working with
(large) systems of equations that arise from the spatial approximation of the
problem. Not surprisingly, the co-existence of a _hierarchy_ of models for a
single physical phenomenon is a common trait of most, if not all, scientific
endeavors.

Precisely because several models of different complexity often exist for one
single physical problem, sometimes it proves convenient to combine several of
them to exploit their relative advantages. For example, when studying complex
deformable bodies, it has been proven effective to combine models for beams
and solids [1, 2, 3], since the economy of the beam equations can be exploited
to analyze whole structures, whereas the (more complex) equations of solid
mechanics are used to describe with detail the mechanics of regions with
intricate stress distributions. This same strategy can be found in the
analysis of complex — typically multiscale — problems in, e.g., the study of
ground and subsurface hydraulic flow [4], arteries and the heart [5],
capillary network and major blood vessels [6, 7], etc. In these situations,
always motivated by a reduction of complexity or computational cost, each of
the connected models might be well-known, but linking them poses difficulties.
In particular, the well-posedness of the joint problem is a delicate matter:
even when each individual model is described with well-posed equations, one
still needs to prove that the connection does not spoil this property. Hence,
links between models of different nature are of theoretical as well as
practical interest.

In this article, we analyze the coupled solution of two diffusive models, the
first one defined over a continuum — and thus described by a PDE — and the
second one defined over discrete network elements — and described with an
ordinary differential equation (ODE). A prototypical example of the problem of
interest in this article is thermal equilibrium. Its most general description
in a three-dimensional body employs Poisson’s equation, the paradigmatic
elliptic model. In addition, when a body is slender, its thermal equilibrium
might be described by a second-order ODE. In practical applications, however,
we might be interested in modeling the temperature of a part that is best
described as a conductive body with a _wire_ that connects some regions, a
wire that need not be inside the body. While one could model the ensemble as a
continuum, it proves more convenient — especially when using numerical
discretizations — to use different models for the bulk and the wire. This is a
relevant problem, for example, in the thermal design of printed circuit boards
(PCBs) where the thermal conductivity of the electronic components and their
(thin) connections are very different, as also their geometries. A second
relevant example that involves coupled models of different nature appears when
modeling the diffusion of infectious diseases. While multi-compartmental
network models have been used for centuries [8] and many results have been
obtained (e.g., [9]), they lack spatial resolution. More recent efforts —
especially related to COVID-19 — attempt to use diffusive boundary value
problems instead [10, 11, 12, 13].

Given its practical and theoretical interests, the numerical approximation of
problems that mix models of different dimensionality has been studied before
(e.g., [14, 3, 15, 16, 6]). In particular, some works have studied the
approximation of mixed-dimensional problems where the low-dimension model is
embedded within the high-dimensional one, with coupling _fluxes_ between them
[14, 6]. One of the main difficulties for coupling 3D-1D models is that the
trace operator from the 3D domain to the 1D domain is not well-posed if the
lower dimensional problem is a 1D manifold with a co-dimension larger than
one. Here, however, we are only interested in coupling elements of a discrete
(network) model with a continuum model. In contrast with the references
before, the network is just a mathematical abstraction, where edges represent
connections (i.e., carriers of flux) and are not, strictly speaking, 1D
submanifolds of the continuum. More specifically, the networks that we will
study in this work model diffusive phenomena between disjoint regions of the
continuum domain, connecting the average values of the linked variables
through a discrete diffusive law. Note that our solution approach differs from
the so-called network diffusion models aiming to solve PDEs on discrete
graphs[17, 18].

Given the difficulties in formulating and analyzing general coupled continuum-
discrete models, we restrict our presentation to the study of the problems
that link two dissimilar models, each of them being the solution of a
minimization problem. The motivation for this choice is double: first, many
diffusive problems of interest have this form and, second, they are easily
amenable to analysis. In particular, if two different models are of this type,
their link might be conveniently tackled using the classical method of
Lagrange multipliers that results in a _saddle point_ problem. The
mathematical aspects of saddle point problems and their approximations with
Galerkin-type methods [19, 20, 21] are well-understood both in
$\mathbb{R}^{n}$ [22] as well as in Hilbert spaces. The current interest
remains, thus, in formulating new links for different models and proving that
the resulting coupled problem as well as its discretization are stable.

In the current work, we study linked formulations of bulk and one-dimensional
diffusive network models, aiming to provide a rigorous footing for the family
of continuum/discrete problems mentioned above. The results presented herein
address first the continuum model, i. e., the well-posedness of the _mixed_
problem that appears when the bulk and one-dimensional discrete type problems
are linked. Once this problem is studied using standard tools, we formulate
mixed finite elements for discrete/continuum diffusion problems and analyze
also their well-posedness. In contrast to elliptic problems, finite element
discretizations of mixed problems do not inherit the well-posedness from the
continuum counterpart and a different stability analysis has to be performed.
In this work, we show the unconditional stability of the continuum/discrete
problems of interest by proving a discrete _inf-sup_ condition. These
theoretical findings are then illustrated with numerical examples. The main
result of this work is thus that linked discrete/continuum finite element
formulations of the coupled diffusive problems considered are unconditionally
stable and convergent. Also, as we will show, the ideas of the proposed
coupling can be used for the solution of diffusive problems with missing data
or partial information. An elegant solution to this seemingly unrelated
problem — closer to data science than diffusion — can be obtained in a
straightforward manner employing the formulation presented in this work.

The article is structured as follows. In Section 2 we describe the
mathematical formulation of the continuum/discrete problems. The theoretical
investigations that prove the well-posedness of the joint problem are
presented in Section 3. Then, in Section 4, we introduce the finite element
discretization of the coupled problem and prove that the well-posedness of the
original problem is inherited by the fully discrete problem. In Section 5, we
demonstrate the applicability of the previously introduced concepts and
methods for several problems that are of practical interest. Finally, Section
6 summarizes the main results of the article.

##  2 Mixed-dimensional Poisson problems

In this section, we define the diffusive problems whose analysis and
approximated solution are the central topic in this article. The first part of
this boundary-value problem describes the stationary solution of a transport
problem in a continuum and is ubiquitous in Mathematical Physics. For
concreteness, in the following, we use the language of thermal analysis
assuming that the unknown field is the temperature in a body. Throughout the
article, all the equations could be reinterpreted in terms, for example, of
mass concentration or electrical charge. The second type of problem refers to
the temperature distribution on a one-dimensional _wire_ whose solution is
given by a second-order _ordinary_ differential equation. In certain
situations, moreover, a closed-form solution to this problem can be found,
yielding a discrete diffusion relation between the temperature at the ends of
the wire. Finally, we describe a joint solution of these two problems when
they are solved simultaneously, that is when we consider a body in thermal
equilibrium with two or more disjoint subsets connected by a thermal wire.

###  2.1 The Poisson problem in a continuum domain

The continuum body where we would like to study its temperature distribution
is a bounded open set $\Omega\subset\mathbb{R}^{d}$ with boundary denoted as
$\partial\Omega$. For simplicity, in what follows, we will restrict to $d=2$,
but no fundamental problem arises in the three-dimensional case. The
temperature in this body is a field $\phi\in H_{0}^{1}(\Omega)$, the Hilbert
space of (Lebesgue) square-integrable functions with square-integrable weak
first derivatives and vanishing trace at the boundary. If $f\in
H^{-1}(\Omega)$ is a known field of heat supply, the stationary Poisson
problem can be written in its standard strong form

$\begin{split}-\kappa\,\Delta\phi&=f\qquad\text{ in }\Omega\ ,\\\
\phi&=0\qquad\text{ on }\partial\Omega\ .\\\ \end{split}$ (1)

Here $\kappa$ is the thermal conductivity of the medium and will be assumed to
be constant, for simplicity, and $\Delta$ is the Laplacian operator. See,
e.g., [23] for a detailed description of this canonical elliptic boundary-
value problem.

With a view to the analysis and discretization of Eq. (1), we rewrite the
Poisson problem in weak form. For that we define $\mathcal{U}\equiv
H_{0}^{1}(\Omega)$ and recall that its norm is given, for any
$\phi\in\mathcal{U}$ by the standard expression

$\|\phi\|_{\mathcal{U}}:=\left(\|\phi\|^{2}_{L^{2}(\Omega)}+\ell^{2}\|\nabla\phi\|^{2}_{L^{2}(\Omega)}\right)^{1/2},$
(2)

where $\ell$ is the characteristic length of the problem (for example, the
diameter of $\Omega$) and $\nabla$ denotes the gradient operator. Then, the
weak form of Poisson’s problem consists in finding $\phi\in\mathcal{U}$ such
that

$a_{\Omega}(\phi,\psi)=f_{\Omega}(\psi)$ (3)

for all $\psi\in\mathcal{U}$, where

$a_{\Omega}(\phi,\psi):=\int_{\Omega}\nabla\phi\cdot\kappa\nabla\psi\;\mathrm{d}V,\qquad
f_{\Omega}(\psi):=\int_{\Omega}f\cdot\psi\;\mathrm{d}A,$ (4)

are, respectively, a bilinear and a linear form on $\mathcal{U}$.

###  2.2 The Poisson problem on a segment

We describe next the stationary heat problem for a one-dimensional body, which
is referred to throughout as a _wire_. The governing equations of this problem
can be derived from the statement (1) of the Poisson problem, simply by
assuming that the body $\Omega$ has a prismatic shape and the temperature
field is constant in all the points of a cross section. Details of this
projection are omitted and the final form of the equation is given.

Consider a wire of length $L$ with temperature $\theta:[0,L]\to\mathbb{R}$.
When the wire is in thermal equilibrium, the temperature must satisfy

$\bar{\kappa}\dfrac{\mathrm{d}^{2}\theta}{\mathrm{d}x^{2}}=r,$ (5)

where $r$ is the heat supply per unit length and $\bar{\kappa}$ is the cross-
sectional conductivity. The problem, however, is not well-posed unless we
append suitable boundary conditions. Rather, the temperature is only defined
modulo an affine function.

As before, we are interested in the weak formulation of this problem.
Postponing for the moment the issue of the uniqueness in the solution, we can
introduce the function space $\mathcal{W}=H^{1}(0,L)$, analogous to the
solution space for the body but without the trace constraint, and look for
solutions $\theta\in\mathcal{W}$ such that

$a_{L}(\theta,\eta)=f_{L}(\eta),$ (6)

for all $\eta\in\mathcal{W}$, where

$a_{L}(\theta,\eta):=\int_{0}^{L}\bar{\kappa}\;\theta^{\prime}\eta^{\prime}\;\mathrm{d}x,\qquad
f_{L}(\eta):=\int_{0}^{L}r\;\eta\;\mathrm{d}x.$ (7)

For future reference, we recall that for every function
$\theta\in\mathcal{W}$, its norm is

$\|\theta\|_{\mathcal{W}}:=\left(\|\theta\|^{2}_{L^{2}(0,L)}+L^{2}\,\|\theta^{\prime}\|^{2}_{L^{2}(0,L))}\right)^{1/2}.$
(8)

Let us note that the actual geometry of the wire plays no role in the
equations above, except for its length. As explained in Sec. 1, Eq. (6) merely
defines the behavior of a network edge, an abstract entity that connects
temperature of two points through a diffusive equation. Moreover, if the heat
supply $r$ is identically zero and the diffusivity is constant, the solution
to Eq. (5) is an affine function. In this case, for all practical purposes,
the wire will just establish a _discrete_ diffusive relation between the
temperature at its two end points.

###  2.3 Linked formulation

We study next a joint problem consisting of a body $\Omega$ with a temperature
field $\phi$ satisfying the problem (3), where additionally we identify two
disjoint non-empty regions $B_{1},B_{2}\subsetneq\Omega$ whose temperatures
are connected by means of a thermally conductive wire. This wire is such that
the temperature at each of its ends coincides with the _mean_ temperature of
the regions $B_{1},B_{2}$, respectively. See Fig. 1 for an illustration of the
linked problem. We note that, alternatively, we could have considered two
disjoint conductive bodies $\Omega_{1}$ and $\Omega_{2}$ that are in thermal
equilibrium while a wire connects regions $B_{\alpha}\subset\Omega_{\alpha}$,
with $\alpha=1,2$. The analysis of both of the problems described is analogous
and, for conciseness, we focus on the former.

Figure 1: Sketch of linked solid-wire problem setting.

Let us insist, once again, that the purpose of this coupled model is not the
actual representation of a true wire that connects two points in the body
(say, the centers of $B_{1}$ and $B_{2}$ in Fig. 1). Rather, the wire connects
regions $B_{1}$ and $B_{2}$ through a discrete diffusive equation. More
precisely, let $\theta_{1},\theta_{2}$ denote the mean temperature in the body
regions $B_{1},B_{2}$, that is,

$\theta_{\alpha}:=\frac{1}{|B_{\alpha}|}\int_{B_{\alpha}}\phi\;\mathrm{d}V,$
(9)

with $\alpha=1,2$, and $|B_{\alpha}|$ denoting the non-zero measure of the
set. Then, the problem governing the temperature field on the wire can be
fully described by the boundary value problem

$\displaystyle\bar{\kappa}\dfrac{d^{2}\theta}{dx^{2}}$
$\displaystyle=0\qquad\text{ on }(0,L),$ (10a) $\displaystyle\theta(0)$
$\displaystyle=\theta_{1}\ ,$ (10b) $\displaystyle\theta(L)$
$\displaystyle=\theta_{2}\ .$ (10c)

If the scalars $\theta_{1},\theta_{2}$ were known, the problem (10) would be
standard and its well-posedness would require no further analysis. However,
here we are interested in the situation where the values
$\theta_{1},\theta_{2}$ are not known _a priori_ but rather, obtained through
the averages (9) of the solution to problem (3).

It remains, thus, to formulate the coupled problem that includes the thermal
equilibrium of the body and wire, as well as the link conditions. We will use
Lagrange multipliers $(\lambda_{1},\lambda_{2})=:\lambda\in
Q\equiv\mathbb{R}^{2}$ to enforce the two constraints (10b) and (10c), and we
will use the notation $\|\cdot\|_{Q}$ to indicate the Euclidean norm in
$\mathbb{R}^{2}$. For convenience, we introduce the product space
$\mathcal{V}=\mathcal{U}\times\mathcal{W}$ with norm

$\|(\phi,\theta)\|_{\mathcal{V}}:=\left(\|\phi\|^{2}_{\mathcal{U}}+\ell^{d-1}\,\|\theta\|^{2}_{\mathcal{W}}\right)^{\frac{1}{2}},$
(11)

for all $(\phi,\theta)\in\mathcal{V}$. On this space we can define the
bilinear form $a(\cdot,\cdot):\mathcal{V}\times\mathcal{V}\to\mathbb{R}$ and
linear form $f:\mathcal{V}\to\mathbb{R}$ as

$\displaystyle a(\phi,\theta;\psi,\delta)$
$\displaystyle:=a_{\Omega}(\phi,\psi)+a_{L}(\theta,\delta)\ ,$ (12)
$\displaystyle f(\psi,\delta)$
$\displaystyle:=f_{\Omega}(\psi)+f_{L}(\delta)\in\mathcal{V}.$

for all $(\phi,\theta)$ and $(\psi,\delta)$ in $\mathcal{V}$.

The joint equilibrium of the solid and the wire then results from the saddle
point of the Lagrangian $\mathcal{L}:\mathcal{V}\times Q\to\mathbb{R}$:

$\begin{split}\mathcal{L}(\phi,\theta,\lambda_{1},\lambda_{2}):=&\frac{1}{2}a(\phi,\theta;\phi,\theta)-f(\phi,\theta)\\\
&+\lambda_{1}\left(\theta(0)-\frac{1}{|B_{1}|}\int_{B_{1}}\phi\;\mathrm{d}V\right)+\lambda_{2}\left(\theta(L)-\frac{1}{|B_{2}|}\int_{B_{2}}\phi\;\mathrm{d}V\right).\end{split}$
(13)

Hence, we aim to solve the following problem

$(\phi,\theta,\lambda_{1},\lambda_{2})^{\star}=\arg\inf_{\phi,\theta}\sup_{\lambda}\mathcal{L}(\phi,\theta,\lambda_{1},\lambda_{2})$
(14)

where the Lagrangian as in Eq. 13. The optimality conditions of the functional
$\mathcal{L}$ are satisfied by the functions
$(\phi,\theta,\lambda_{1},\lambda_{2})\in\mathcal{V}\times Q$ such that

$\displaystyle
a(\phi,\theta;\psi,\delta)+b(\psi,\delta;\lambda_{1},\lambda_{2})$
$\displaystyle=f(\psi,\delta)\ ,$ (15) $\displaystyle
b(\phi,\theta;\mathit{\Gamma}_{1},\mathit{\Gamma}_{2})$ $\displaystyle=0,$
(16)

for all
$(\psi,\delta,\mathit{\Gamma}_{1},\mathit{\Gamma}_{2})\in\mathcal{V}\times Q$
with

$b(\phi,\theta;\mathit{\Gamma}_{1},\mathit{\Gamma}_{2}):=\mathit{\Gamma}_{1}\left(\theta(0)-\frac{1}{|B_{1}|}\int_{B_{1}}\phi\,\;\mathrm{d}V\right)+\mathit{\Gamma}_{2}\left(\theta(L)-\frac{1}{|B_{2}|}\int_{B_{2}}\phi\,\;\mathrm{d}V\right).$
(17)

The bilinear form $a(\cdot,\cdot)$ defines a linear continuum operator
$A:\mathcal{V}\to\mathcal{V}^{\prime}$ by the relation

$\langle Au,v\rangle_{\mathcal{V}^{\prime}\times\mathcal{V}}=a(u,v),$ (18)

for all $u=(\phi,\theta)\in\mathcal{V},v=(\psi,\delta)\in\mathcal{V}$. Also,
the bilinear form $b(\cdot,\cdot)$ on $\mathcal{V}\times Q$ defines a linear
operator $B:\mathcal{V}\to Q^{\prime}$ with transpose
$B^{T}:Q\to\mathcal{V}^{\prime}$ by

$\langle Bv,\mathit{\Gamma}\rangle_{Q^{\prime}\times Q}=\langle
v,B^{T}\mathit{\Gamma}\rangle_{\mathcal{V}\times\mathcal{V}^{\prime}}=b(v,\mathit{\Gamma}),$
(19)

for all $v\in\mathcal{V},\mathit{\Gamma}\in Q$. Employing these definitions,
Eq. (16) can be alternatively rewritten as:

$\displaystyle Au+B^{T}\lambda$ $\displaystyle=f\qquad\text{ in
}\mathcal{V}^{\prime}$ (20) $\displaystyle Bu$ $\displaystyle=0\qquad\text{ in
}Q^{\prime}.$

This last expression is the standard form of a _mixed problem_ [20]. The
solvability of this problem depends on conditions over the bilinear forms
$a(\cdot,\cdot)$ and $b(\cdot,\cdot)$ as well as properties of the spaces
where they are defined on. In particular, often the properties of the linear
operator $B$ are delicate to ascertain.

######  Remarks 1.

Two modifications of problem (16) are interesting in their own right:

  1. 1.

One could consider that the one-dimensional diffusive mechanism connects,
instead of disjoint subsets $B_{1},B_{2}\subset\Omega$, the boundary of two
different bodies, or parts of them. The description of this modified problem
and the functional setting are slightly different than the one presented up to
here, since the new problem will be formulated in terms of the traces of
functions.

  2. 2.

Second, we could consider a situation where there is no connecting one-
dimensional diffusive wire between regions of the body but only that the
_average_ temperature on a measurable set $S\subset\Omega$ is known to have a
fixed value $\bar{\theta}$. In this simple case we will be left with Poisson’s
problem with a constraint. The problem will still be of the form (16), more
precisely,

$\displaystyle a_{\Omega}(\phi,\psi)+\bar{b}(\psi,\gamma)$
$\displaystyle=f_{\Omega}(\psi)\ ,$ (21) $\displaystyle\bar{b}(\phi,\eta)$
$\displaystyle=0,$

with $\bar{b}(\cdot,\cdot)$ on $\mathcal{V}\times\mathbb{R}$ defined as

$\bar{b}(\phi,\eta):=\eta\left(\bar{\theta}-\frac{1}{|S|}\int_{S}\phi\,\mathrm{d}V\right).$
(22)

##  3 Analysis

In this section, we study the well-posedness of problem (16). According to the
standard theory of mixed problems [19, 20, 21], we need to show that the
bilinear forms $a(\cdot,\cdot)$ and $b(\cdot,\cdot)$ are continuum, that
$a(\cdot,\cdot)$ is elliptic on the kernel of the operator $B$, and that a
certain _inf-sup_ condition, to be defined later, also holds for
$b(\cdot,\cdot)$.

In the following, for simplicity, let us assume that the conductivities
$\kappa$ and $\bar{\kappa}$ are constant and positive. The first step of the
analysis is to study the continuity of all the linear and bilinear forms
appearing in the problem statement (16). This is trivial and we summarize all
the results, without proof, in the following theorem:

######  Theorem 3.1.

The bilinear forms $a_{\Omega}(\cdot,\cdot)$, $a_{L}(\cdot,\cdot)$, and
$a(\cdot,\cdot)$ are continuous in their corresponding spaces of definition,
i.e.,

$\displaystyle|a_{\Omega}(\phi,\psi)|$ $\displaystyle\leq
c_{\Omega}\;\|\phi\|_{\mathcal{U}}\;\|\psi\|_{\mathcal{U}}\ ,$ (23)
$\displaystyle|a_{L}(\theta,\delta)|$ $\displaystyle\leq
c_{L}\;\|\theta\|_{\mathcal{W}}\;\|\delta\|_{\mathcal{W}}\ ,$
$\displaystyle|a(\phi,\theta;\psi,\delta)|$ $\displaystyle\leq
c\;\|(\psi,\delta)\|_{\mathcal{V}}\ ,$

for all $\phi,\psi\in\mathcal{U}$, $\theta,\delta\in\mathcal{W}$ and some
generic positive constants $c_{\Omega},c_{L}$ and $c$.

Next, we show the continuity of $b(\cdot,\cdot)$. Since this bilinear form is
non-standard, we provide the full proof of the result.

######  Theorem 3.2.

The bilinear form $b$ is continuous on $(\mathcal{V}\times Q)$, i.e.

$|b(\phi,\theta;\lambda)|\leq
c\;\|(\phi,\theta)\|_{\mathcal{V}}\;\|\lambda\|_{Q}$ (24)

for all $(\phi,\theta;\lambda)\in\mathcal{V}\times Q$ and some constant $c>0$.

###### Proof.

Since $H^{1}(0,L)\hookrightarrow C^{0}[0,L]$, we can use the mean value
theorem to determine that there exists $m\in[0,L]$ such that

$\theta(m)=\frac{1}{L}\int_{0}^{L}\theta(x)\;\mathrm{d}x\ .$ (25)

Then, by the fundamental theorem of calculus

$\theta(L)=\theta(m)+\int_{m}^{L}\theta^{\prime}(x)\;\mathrm{d}x.$ (26)

Hence,

$\displaystyle|\theta(L)|$
$\displaystyle\leq\frac{1}{L}\int_{0}^{L}|\theta(x)|\;\mathrm{d}x+\int_{m}^{L}|\theta^{\prime}(x)|\;\mathrm{d}x$
(27) $\displaystyle\leq
L^{-1/2}\,\|\theta\|_{L^{2}(0,L)}+{L}^{1/2}\,\|\theta^{\prime}\|_{L^{2}(0,L)}$
$\displaystyle={L}^{-1/2}\,\|\theta\|_{H^{1}(0,L)}.$

Similarly, the bound $|\theta(0)|\leq{L}^{-1/2}\,\|\theta\|_{H^{1}(0,L)}$ also
holds. Using these results, we proceed to bound from above the bilinear form
$b(\cdot,\cdot)$:

$\displaystyle|b(\phi,\theta;\lambda)|$
$\displaystyle=\left|\lambda_{1}\left(\frac{1}{|B_{1}|}\int_{B_{1}}\phi\;\mathrm{d}V-\theta(0)\right)+\lambda_{2}\left(\frac{1}{|B_{2}|}\int_{B_{2}}\phi\;\mathrm{d}V-\theta(L)\right)\right|$
$\displaystyle\leq|\lambda_{1}|\left|\frac{1}{|B_{1}|}\int_{B_{1}}\phi\;\mathrm{d}V-\theta(0)\right|+|\lambda_{2}|\left|\frac{1}{|B_{2}|}\int_{B_{2}}\phi\;\mathrm{d}V-\theta(L)\right|$
$\displaystyle\leq|\lambda_{1}|\left(\frac{1}{|B_{1}|}\int_{B_{1}}|\phi|\;\mathrm{d}V+|\theta(0)|\right)+|\lambda_{2}|\left(\frac{1}{|B_{2}|}\int_{B_{2}}|\phi|\;\mathrm{d}V+|\theta(L)|\right)$
$\displaystyle\leq
C\,\|\lambda\|_{2}\left(\frac{1}{|\Omega|}\int_{\Omega}|\phi|\;\mathrm{d}V+L^{-1/2}\,\|\theta\|_{H^{1}(0,L)}\right)$
$\displaystyle\leq
C\,\|\lambda\|_{2}\left(\ell^{-d/2}\|\phi\|_{H_{0}^{1}(\Omega)}+L^{-1/2}\|\theta\|_{H^{1}(0,L)}\right)$
$\displaystyle\leq
C\;\ell^{-d/2}\|\lambda\|_{Q}\;\|(\phi,\theta)\|_{\mathcal{V}},$

where, throughout the proof, $C$ denotes a constant whose value might change
from one step to another and $d$ denotes the dimension of the characteristic
length corresponding to $\Omega$. ∎

Next, we have to ensure ellipticity on $\ker B$, the kernel of the operator
$B$. By definition, this set is

$\ker
B:=\\{(\phi,\theta)\in\mathcal{V}:\;b(\phi,\theta;\lambda)=0,\quad\forall\lambda\in
Q\\}.$ (28)

Elements $(\phi,\theta)\in\mathcal{V}$ in $\ker B$ verify

$\lambda_{1}\left(\dfrac{1}{|B_{1}|}\int_{B_{1}}\phi\;\mathrm{d}V-\theta(0)\right)+\lambda_{2}\left(\dfrac{1}{|B_{2}|}\int_{B_{2}}\phi\;\mathrm{d}V-\theta(L)\right)=0.$
(29)

for all pairs $(\lambda_{1},\lambda_{2})\in Q$. Since the two Lagrange
multipliers are independent, it must hold that

$\ker
B=\left\\{(\phi,\theta)\in\mathcal{V}:\dfrac{1}{|B_{1}|}\int_{B_{1}}\phi\;\mathrm{d}V=\theta(0)\
\text{ and }\
\dfrac{1}{|B_{2}|}\int_{B_{2}}\phi\;\mathrm{d}V=\theta(L)\right\\}.$ (30)

As advanced above, for the well-posedness of the mixed problem, it suffices
that the bilinear form $a(\cdot,\cdot)$ be elliptic on $\ker
B\subset\mathcal{V}$, as shown in the following theorem:

######  Theorem 3.3.

The bilinear form $a(\cdot,\cdot)$ is elliptic on $\ker B$, i.e., there exists
a constant $\bar{\alpha}>0$ such that

$a(\phi,\theta;\phi,\theta)\geq\bar{\alpha}\;\|(\phi,\theta)\|^{2}_{\mathcal{V}},$
(31)

for all $(\phi,\theta)\in\ker B$.

###### Proof.

To start the proof, we need two preliminary results. First, using the weighted
Young inequality we note that

$\displaystyle\|\theta-$
$\displaystyle\theta(0)\|^{2}_{L^{2}(0,L)}+\|\theta-\theta(L)\|^{2}_{L^{2}(0,L)}$
(32) $\displaystyle=$ $\displaystyle
2\,\|\theta\|^{2}_{L^{2}(0,L)}+L\,|\theta(0)|^{2}+L\,|\theta(L)|^{2}-2\int_{0}^{L}\theta(0)\,\theta\;\mathrm{d}x-2\int_{0}^{L}\theta(L)\,\theta\;\mathrm{d}x$
$\displaystyle\geq$ $\displaystyle
2\,\|\theta\|^{2}_{L^{2}(0,L)}+L\,|\theta(0)|^{2}+L\,|\theta(L)|^{2}$
$\displaystyle-\epsilon_{0}^{2}\,\|\theta\|^{2}_{L^{2}(0,L)}-\frac{L}{\epsilon_{0}^{2}}|\theta(0)|^{2}-\epsilon_{L}^{2}\,\|\theta\|^{2}_{L^{2}(0,L)}-\frac{L}{\epsilon_{L}^{2}}|\theta(L)|^{2}$
$\displaystyle\geq$
$\displaystyle\left(2-\epsilon_{0}^{2}-\epsilon_{L}^{2}\right)\|\theta\|^{2}_{L^{2}(0,L)}+L\left(1-\frac{1}{\epsilon_{0}^{2}}\right)|\theta(0)|^{2}+L\left(1-\frac{1}{\epsilon_{L}^{2}}\right)|\theta(L)|^{2},$

for arbitrary scalars $\epsilon_{0},\epsilon_{L}>0$. In the kernel of $B$,
moreover, we have that

$\displaystyle|\theta(0)|^{2}$
$\displaystyle=\frac{1}{|B_{1}|}\left|\int_{B_{1}}\phi\;\mathrm{d}V\right|^{2}\leq
C\,\ell^{-d}\|\phi\|^{2}_{H^{1}_{0}(\Omega)},$ (33)
$\displaystyle|\theta(L)|^{2}$
$\displaystyle=\frac{1}{|B_{2}|}\left|\int_{B_{2}}\phi\;\mathrm{d}V\right|^{2}\leq
C\,\ell^{-d}\|\phi\|^{2}_{H^{1}_{0}(\Omega)}\ .$

Assuming that the following conditions hold

$1-\frac{1}{\epsilon_{0}^{2}}\leq 0\ ,\qquad 1-\frac{1}{\epsilon_{L}^{2}}\leq
0\ ,\qquad$ (34)

then, combining Eqs. 32, LABEL: and 33, the following bound is obtained

$\displaystyle\|\theta-$
$\displaystyle\theta(0)\|^{2}_{L^{2}(0,L)}+\|\theta-\theta(L)\|^{2}_{L^{2}(0,L)}$
(35) $\displaystyle\geq$
$\displaystyle\left(2-\epsilon_{0}^{2}-\epsilon_{L}^{2}\right)\|\theta\|^{2}_{L^{2}(0,L)}+C\,\ell^{1-d}\left(2-\frac{1}{\epsilon_{0}^{2}}-\frac{1}{\epsilon_{L}^{2}}\right)\|\phi\|^{2}_{H^{1}_{0}(\Omega)}\
.$

Also, as in the proof of Theorem 3.2, we use the fundamental theorem of
calculus to bound

$\displaystyle\|\theta^{\prime}\|^{2}_{L^{2}(0,L)}$ $\displaystyle\geq
CL^{-2}\|\theta-\theta(0)\|^{2}_{L^{2}(0,L)}\ ,$ (36)
$\displaystyle\|\theta^{\prime}\|^{2}_{L^{2}(0,L)}$ $\displaystyle\geq
CL^{-2}\|\theta-\theta(L)\|^{2}_{L^{2}(0,L)}\ .$

To finally prove the ellipticity of the bilinear form $a(\cdot,\cdot)$ we use
the definition (12), the bounds (35), (36) and the Poincaré inequality from
where it follows that

$\displaystyle a(\phi,\theta;\phi,\theta)=$
$\displaystyle\int_{\Omega}\kappa\;|\nabla\phi|^{2}\;\mathrm{d}V+\int_{0}^{L}\bar{\kappa}\;|\theta^{\prime}|^{2}\;\mathrm{d}x$
$\displaystyle\geq$ $\displaystyle
C\,\kappa\,\ell^{-2}\|\phi\|^{2}_{H^{1}_{0}(\Omega)}+\frac{\bar{\kappa}}{2}\|\theta^{\prime}\|^{2}_{L^{2}(0,L)}+\frac{\bar{\kappa}}{2}\|\theta^{\prime}\|^{2}_{L^{2}(0,L)}$
$\displaystyle\geq$ $\displaystyle
C\,\kappa\,\ell^{-2}\|\phi\|^{2}_{H^{1}_{0}(\Omega)}+C\,\kappa\,\ell^{d-1}\|\theta^{\prime}\|^{2}_{L^{2}(0,L)}$
$\displaystyle+C\,\kappa\,\ell^{d-3}\left(\|\theta-\theta(0)\|^{2}_{L^{2}(0,L)}+\|\theta-\theta(L)\|^{2}_{L^{2}(0,L)}\right)$
$\displaystyle\geq$ $\displaystyle
C\,\kappa\,\ell^{-2}\|\phi\|^{2}_{H^{1}_{0}(\Omega)}+C\,\kappa\,\ell^{d-1}\|\theta^{\prime}\|^{2}_{L^{2}(0,L)}$
$\displaystyle+C\,\kappa\,\ell^{d-3}\left(2-\epsilon_{0}^{2}-\epsilon_{L}^{2}\right)\|\theta\|^{2}_{L^{2}(0,L)}$
$\displaystyle+C\,\kappa\,\ell^{-2}\left(2-\frac{1}{\epsilon_{0}^{2}}-\frac{1}{\epsilon_{L}^{2}}\right)\|\phi\|^{2}_{H^{1}_{0}(\Omega)},$

where, as usual, $C$ denotes a generic constant whose value might change in
each inequality. Finally, we rewrite this bound as

$\displaystyle a(\phi,\theta;\phi,\theta)\geq$ $\displaystyle
C\,\kappa\,\ell^{-2}\left(3-\frac{1}{\epsilon_{0}^{2}}-\frac{1}{\epsilon_{L}^{2}}\right)\|\phi\|^{2}_{H^{1}_{0}(\Omega)}$
$\displaystyle+C\,\kappa\,\ell^{d-3}\left((2-\epsilon_{0}^{2}-\epsilon_{L}^{2})\|\theta\|^{2}_{L^{2}(0,L)}+L^{2}\|\theta^{\prime}\|^{2}_{L^{2}(0,L)}\right)$
$\displaystyle\geq$ $\displaystyle
C\,\kappa\,\ell^{-2}\left(\|\phi\|^{2}_{\mathcal{U}}+\ell^{d-1}\|\theta\|^{2}_{\mathcal{W}}\right)$
$\displaystyle=$ $\displaystyle
C\,\kappa\,\ell^{-2}\|(\phi,\theta)\|^{2}_{\mathcal{V}}\ ,$

as long as there exist $\epsilon_{0}^{2},\epsilon_{L}^{2}\geq 0$ that verify
conditions (34) as well as

$3-\frac{1}{\epsilon_{0}^{2}}-\frac{1}{\epsilon_{L}^{2}}>0\ ,\qquad
2-\epsilon_{0}^{2}-\epsilon_{L}^{2}>0\ .$

These three conditions are satisfied, for example, by
$\epsilon_{0}^{2}=\epsilon_{L}^{2}=4/5$, completing the proof. ∎

Let us note that the bilinear form $a(\cdot,\cdot)$ is not elliptic in the
whole space $\mathcal{V}$, because the bilinear form of the one-dimensional
conductor, namely $a_{L}(\cdot,\cdot)$ is not elliptic in the space
$\mathcal{W}$, precisely due to the lack of Dirichlet boundary conditions in
problem (5).

Finally, to ensure that the mixed problem is well-posed, we also have to show
that the following _inf-sup_ condition holds for $b(\cdot,\cdot)$.

######  Theorem 3.4.

Inf-sup condition. There exists a constant $\beta>0$ such that

$\sup_{(\phi,\theta)\in\mathcal{V}}\frac{b(\phi,\theta;\lambda)}{\|(\phi,\theta)\|_{\mathcal{V}}}\geq\beta\|\lambda\|_{Q\backslash\ker
B^{T}}.$ (37)

###### Proof.

Since $B$ is an operator from $\mathcal{V}$ to $Q^{\prime}$, a finite-
dimensional space, its range is closed (Theorem 1.1. in [20]). Hence, the inf-
sup condition holds. ∎

With Theorems 3.1, LABEL:, 3.2, LABEL:, 3.3, LABEL: and 3.4, it follows that
problem (16) is well-posed. Let us remark that, typically, in other
constrained minimization problems, the proof of the _inf-sup_ condition might
be quite involved. In the problem discussed in this article, however, this
proof is trivial.

##  4 Approximation of the problem

To solve the saddle-point problem (16) we use a mixed finite element
discretization. We now briefly recall the details of this method as it applies
to the project at hand.

The first step in the finite element approximation of problem (16) is the
definition of a finite-dimensional subspace, $\mathcal{V}_{h}\subset V$ of the
infinite-dimensional solution space. Moreover, we assume that all the finite
element functions in $\mathcal{V}_{h}$ are linear combinations of piecewise
polynomials defined in their corresponding domains, namely, $\Omega$ or
$[0,L]$. Since the space of the Lagrange multipliers is already a space of
dimension 2 we do not need to introduce a subspace for it and we simply define
$Q_{h}\equiv Q$.

The (mixed) Galerkin finite element method consists in finding
$u_{h}=(\phi_{h},\theta_{h})\in\mathcal{V}_{h},\;\lambda_{h}\in Q_{h}$ such
that

$\displaystyle a(u_{h},v_{h})+b(v_{h},\lambda_{h})$ $\displaystyle=f(v_{h}),$
(38) $\displaystyle b(u_{h},\mathit{\Gamma}_{h})$ $\displaystyle=0,$

holds for all $v_{h}=(\psi_{h},\delta_{h})\in\mathcal{V}_{h}$ and
$\mathit{\Gamma}_{h}\in Q_{h}$. In contrast with the finite element
approximations of elliptic boundary value problems, the well-posedness of
mixed methods such as (38) does not follow from the well-posedness of the
corresponding continuous problem. Instead, a new analysis has to be carried
out and, in particular, a discrete _inf-sup_ condition needs to be proven.
This is typically the most difficult ingredient of this analysis but, as we
will show below, it is not the case for the problem at hand given the finite
dimension of the space $Q$.

###  4.1 Analysis

The well-posedness of the discrete problem (38) can be established by
following similar steps as in the analysis of the continuous problem and
presented in Section 3. For that, we start by introducing $B_{h}$, the
discrete counterpart of the operator $B$ used in Section 3, and now defined by
the relation

$\langle B_{h}v_{h},\mathit{\Gamma}_{h}\rangle_{Q_{h}^{\prime}\times
Q_{h}}=\langle
v_{h},B^{T}_{h}\mathit{\Gamma}_{h}\rangle_{\mathcal{V}\times\mathcal{V}^{\prime}}=b(v_{h},\mathit{\Gamma}_{h}),$
(39)

for all $v_{h}\in\mathcal{V}_{h}$ and $\mathit{\Gamma}_{h}\in Q_{h}$. Since
$B_{h}$ is surjective and $\ker\mathcal{B}_{h}^{T}\subset\ker\mathcal{B}^{T}$,
$\ker B_{h}^{T}=\\{0\\}$ and $B_{h}=B{\big{|}}_{V_{h}}$. So, we end up in the
special case of $\ker B_{h}\subset\ker B$. Hence, the ellipticity of
$a(\cdot,\cdot)$ on $\ker B_{h}$ follows from the ellipticity of
$a(\cdot,\cdot)$ on $\ker B$. As advanced, the key condition for the well-
posedness of (38) is thus the _discrete inf-sup_ condition. However, given
that the range of $B_{h}$ is finite dimensional, it is closed and thus the
_discrete inf-sup_ condition holds. From Theorem 2.1. in [20] we obtain the
following estimate

$\|u-u_{h}\|_{\mathcal{V}}+\|\lambda-\lambda_{h}\|_{Q\backslash\ker{B^{T}}}\leq
c\left(\inf_{v_{h}\in\mathcal{V}_{h}}\|u-v_{h}\|_{\mathcal{V}}+\inf_{q_{h}\in
Q_{h}}\|\lambda-\mathit{\Gamma}_{h}\|_{Q}\right),$ (40)

where $c$ is a constant that depends on $a(\cdot,\cdot)$, $b(\cdot,\cdot)$,
$\bar{\alpha}$ as in Theorem 3.3 and $\beta$ as in Theorem 3.4. Details can be
found in [20].

##  5 Numerical examples

In this section, we use first the mixed method (38) for two examples chosen to
illustrate its ability to link regions of thermally conductive solids. Then,
and to complete the section, we show that the ideas of the mixed network-
continuum formulations can be exploited to solve one particular type of
inference problems in diffusive situations.

###  5.1 Heat flux between dissimilar regions

In this example, we choose a square solid with dimensions $2\times 2$ where
two subregions, the first with the shape of a circle, and the second an
ellipse, are linked through a conductive wire. In a Cartesian coordinate
system located at the center of the square and with axes parallel to its
edges, the circular region has center $(x,y)=(-0.625,-0.5625)$ and radius
$0.5$. The elliptical region is centered at $(x,y)=(0.375,0.4375)$, and has
horizontal and vertical semi-axes of lengths $0.175$ and $0.3$, respectively
(see Fig. 2). The setting is very similar to the one we employed in Section
2.3 to describe the mixed-dimensional boundary value problem. See Fig. 1 for
comparison.

The block has conductivity $\kappa=100$ and the wire has a high conductivity
of value $\bar{\kappa}=10,000$. The temperature is constrained at the top and
bottom edges to values of $500$ and $300$, respectively. If the problem had no
thermal link, the thermal field would be linear with constant vertical heat
flux. However, due to the presence of the linking wire, which is selected with
a high conductivity, the thermal field is distorted. The elliptical region
–closer to the hotter top– serves as a heat source for the circular region
–this one closer to the cooler bottom edge. As a result, the temperature
fields in the two connected regions are close to 400.

In principle, the region and its center can lay anywhere in the domain. The
centers of these regions determine the amount of temperature going into and
coming out of the wire calculated as the average temperature in the
corresponding region. To study the convergence of the formulation, we obtain
six finite element solutions for this problem, with increasing resolution (see
Fig. 2). In each of the solutions, the linked regions $B_{1},B_{2}$ consist of
those elements of the mesh whose centers fall inside the circle and ellipse,
respectively. As can be observed in Fig. 2, coarse meshes do not represent
accurately the linked regions (see the top two figures in Fig. 2). However,
when the mesh is refined, both the circle as well as the ellipse are
accurately approximated (see, again, Fig. 2). We note that, alternatively, we
could have chosen to employ a mesh that was adapted to the linked regions, but
the converged solutions would not differ.

(a) $8^{2}$ mesh.

(b) $16^{2}$ mesh.

(c) $32^{2}$ mesh.

(d) $64^{2}$ mesh.

(e) $128^{2}$ mesh.

(f) $256^{2}$ mesh.

Figure 2: Thermal link connecting dissimilar regions of circle and ellipse
inside a block for six meshes with $8^{2}$, $16^{2}$, $32^{2}$, $64^{2}$,
$128^{2}$, $256^{2}$ elements (top left to bottom right).

###  5.2 Convergence of increasingly finer connected regions

This second example examines the classical Fourier problem and its relation
with diffusive networks, of the type employed in this work, with the continuum
notion of flux. For that, we will consider a square domain of unit side, a
material with conductivity $\kappa=100$, and temperature boundary conditions
described by parabolas at the four edges, all with their maximum values at the
center of the edges and zero at the vertices. The maximum temperatures at each
edge are $100$ (bottom), $200$ (right), $300$ (top), and $400$ (left).

This problem can be solved with a classical numerical method, such as the
finite element or finite volume method (see Fig. 4(a) for the finite element
solution with a mesh of $512^{2}$ bilinear elements). However, here we solve
the problem by discretizing the domain into disjoint and _independent_ square
regions. These regions are connected with their neighbors through wires with
conductivity $\kappa\cdot\Delta y$ (horizontal wires) and $\kappa\cdot\Delta
x$ (vertical wires), where $\Delta x$ and $\Delta y$ are the lengths of the
corresponding region in $x$ and $y$ direction, respectively. Then, this linked
problem is solved using the formulation described in Section 4. Since we use
only one finite element for each region, the method explained above is
equivalent to a finite volume method where the fluxes are obtained by solving
the boundary value problem of the wires.

Fig. 3 shows the solutions obtained with an increasing number of independent
regions. In addition to the thermal field, the edges of the diffusive network
are also drawn as well as their temperature field, using the same color scale
as in the continuum regions. In the coarser solutions, one can clearly see the
discontinuity of the thermal field across the boundaries of the regions since,
as explained, they are independent meshes and do not share any node. Also, it
is apparent in these coarser solutions that the thermal field in the wires
differs from the continuum field at corresponding points.

Remarkably, as the size of the regions is reduced –and hence the length and
conductivity of the wires– the connected solution (Fig. 3) resembles more
closely the finite element solution (Fig. 4(a)) used as reference. This
apparent convergence is verified by the results depicted in Fig. 4(b). This
figure shows the root-mean-square error of the temperature field of each
linked solution compared with the reference finite element solution and
obtained as

$e=\sqrt{\frac{1}{N}\sum_{i}^{N}(\phi_{i}-\phi_{i}^{FE})^{2}}\ ,$ (41)

where $N$ is the number of regions and $\phi_{i}$ refers to the temperature at
the center of the $i-$th region, or element, respectively.

(a) $2^{2}$ regions

(b) $4^{2}$ regions

(c) $8^{2}$ regions

(d) $16^{2}$ regions

(e) $32^{2}$ regions

(f) $64^{2}$ regions

Figure 3: Temperature field in thermal problem solved with independent regions
linked with wires. The linking wires are depicted as lines.

(a) Finite element solution of the several regions and wires problem with a
quadratic mesh of $512\times 512$ elements.

(b) Root-mean-square error of the regions’ center temperature compared to the
finite element solution. Orange triangle represents linear convergence.

Figure 4: Finite element solution of the several regions and wires problem and
root-mean-square error of the regions’ center temperature compared to the
finite element solution.

###  5.3 Using linked models to infer diffusive solutions

This final example studies the possibility of using the methods introduced in
Section 4 to infer the optimal thermal field in a domain when only partial
information is available. More precisely, we are interested in determining the
most likely temperature field in a domain when there is information about the
temperature on part of the boundary and the _average_ temperature in some
regions of the domain.

When studying a Fourier-type diffusive problem, one needs to know the
Dirichlet and Neumann conditions in all of the boundaries as well as the heat
applied in the interior of the domain, see Eq. 1. The lack of either of these
data renders the problem ill-posed and no solution can be analytically (nor
numerically) obtained. Here we are interested in finding the optimal
temperature field of a continuum domain, $\Omega$, where the Dirichlet
boundary conditions are partially known but also the average temperature in
some regions,
$\\{B_{i}\\}_{i=1}^{N_{regions}}\;\mathrm{with}\;B_{i}\subseteq\Omega$. In
particular, the heat supplied through the Neumann boundary and the interior of
the domain is completely unknown. This is an optimization problem that
searches for the temperature field $\phi\in\mathcal{U}$ that verifies

$\displaystyle\min_{\phi\in\mathcal{U}}$
$\displaystyle\mathcal{E}=\int_{\Omega}\frac{\kappa}{2}\,\|\nabla\phi\|^{2}\mathrm{d}V,$
(42) $\displaystyle\mathrm{s.t.}$
$\displaystyle\frac{1}{|B_{i}|}\int_{B_{i}}\phi\;\mathrm{d}V=\theta_{i},\quad
i\in\\{1,...,N_{regions}\\}$
$\displaystyle\phi=\bar{\phi},\quad\emptyset\neq\partial\Omega_{j}\subseteq\partial\Omega$

The solution to this optimization problem is the thermal field on the solid
$\Omega$ with imposed mean-temperature constraints using _degenerated wires_
that link the assumed temperature with the selected regions using the
constraint given in Eqs. 22 and 22.

To analyze the ability of problem (42) to infer an unknown thermal field, we
consider a $2\times 1$ rectangle with $\kappa=100$. We compute a reference
solution with boundary conditions and thermal loading depicted in Fig. 5(a).
In this figure, the temperature on the left and top edges corresponds to
$\phi=400$ and $\phi=200$, respectively. Also, a wavy heat supply is imposed
in the whole domain of the form

$h(x,y)=\begin{cases}0,&\quad\mathrm{if}\;(x,y)\in\mathcal{C}\\\
\sin\left(\frac{4\pi(x-x_{x})^{2}+(y-y_{c})^{2}}{r^{2}}\right),&\quad\mathrm{otherwise}\end{cases},$
(43)

where $\mathcal{C}$ is the circle with center $(x,y)=(1.5,0.25)$ and radius
$r=0.25$, where the coordinates refer to a Cartesian system located at the
bottom left corner of the rectangle, with the $x,y$ axes parallel to the
horizontal and vertical directions, respectively. Fig. 5(b) shows the solution
of this boundary value problem (in what follows, the reference problem)
obtained with a finite element mesh of size $128\times 64$.

For the current example, we define the potentially linked regions as the ones
obtained with a simple grid of the whole domain of $8\times 8$ equal
rectangles. We only employ as known boundary data the temperature on the left
edge. Then, we analyze a sequence of discrete problems that employ an
increasingly large number of data concerning the average temperature in random
regions. Denoting these regions as $B_{i}$ (see Eq. 42), we proceed as
follows: first, we solve the problem with known Dirichlet data, and a known
average temperature $\theta_{i}$ on a random region $B_{1}$ which we can
obtain from the exact solution; then, we select a second random region
$B_{2}$, obtain its average temperature $\theta_{2}$ from the exact solution
and solve the minimization problem with the two constraints. Proceeding in
this fashion, we obtain the solutions illustrated in Fig. 6. This figure shows
that when the only available data is the average temperature of one region and
the Dirichlet boundary condition, the solution of minimal energy is far from
the reference (see top figures in Fig. 6. As expected, the more data
available, the better the obtained solution. When the average temperatures are
available in all the rectangular regions, the error is minimal and the
solution of the optimization problem is close to the reference solution.
Naturally, the small oscillations in the solution are not captured by the
approximation because the available data do not have enough resolution. One
might expect that in the limit when the regions of known average become very
small and cover the whole domain, the solution of the optimization problem
converges to the exact solution.

Finally, in Fig. 7 the energy error $\varepsilon$ — relative to the energy of
the reference solution —, i.e.,

$\varepsilon=\frac{|\mathcal{E}^{N_{regions}}-\mathcal{E}^{FE}|}{\mathcal{E}^{FE}}$
(44)

for the calculated sequence is illustrated. The energy error depends on the
random sequence of known data, but it should be always monotonically
decreasing. Note that for this sequence, there is almost no improvement in the
solution for the last $33$ added data instances. The average temperature of
the region added in the $32nd$ iteration makes the difference, reducing the
energy error to approximately $20\%$. From that iteration, the solutions of
the constrained optimization problems are very similar. (See Figs. 6(h) and
6(j)).

(a) Heat supply field and boundary conditions.

(b) Reference solution with a fe mesh of size $128\times 64$.

Figure 5: Reference boundary value problem solution of the inference example.

(a) $1$ random region

(b) Solution with $1$ random region

(c) $9$ random regions

(d) Solution with $9$ random regions

(e) $25$ random region

(f) Solution with $25$ random regions

(g) $49$ random regions

(h) Solution with $49$ random regions

(i) All regions

(j) Solution with all regions

Figure 6: Inference problem. A sequence of solutions with a different number
of available region data. Left: gray regions are the ones where the average
temperature is available. Right: solution with the corresponding average
temperatures and only left edge boundary condition known. Figure 7: Energy
error for the inference problem for different numbers of available region
data.

##  6 Conclusions and main results

In this paper, we studied linked continuum/discrete models of diffusion from
the theoretical and numerical points of view. More specifically, we consider
the solution of diffusion boundary value problems where (discrete) diffusion
pathways are introduced between subdomains, more precisely, between the mean
values of the diffusing quantity. Although we have described all the results
and examples in the language of thermal problems, the scope of the article is
more general and applies to any type of linear Fourier-type boundary value
problem.

After introducing the setting for the linked formulations of bulk and the one-
dimensional diffusive models, we proved the well-posedness of the
continuum/discrete linked problem in the corresponding functional spaces by
resorting to the theory of constrained boundary value problems. Next, we
proved that the well-posedness continues to hold for the approximate problem,
i.e. the problem discretized via mixed finite elements. A direct consequence
of the two previous results is the convergence of the fully discretized
equations of the linked model to its exact solution.

The previous theoretical findings have been illustrated with three numerical
examples. They show that the finite element method is always stable and that
the links take care of the diffusive fluxes in a different scale than the
continuum, but that in the limit both discrete and continuum mechanisms are
equivalent. Finally, we showed that the ideas of the method can be used to
carry out inference in the solution of diffusive problems with only partial
information available about the _average_ solution in subsets of the whole
domain.

The results of this article can be employed to link the two different classes
of diffusive models, namely, the PDE-based descriptions of the continuum and
diffusive network-based formulations when the network connects average values
of the diffusive field computed over bounded disjoint regions.

## Acknowledgements

C.S. and I.R. acknowledge the financial support from Madrid’s regional
government through a grant with IMDEA Materials addressing research activities
on SARS-COV 2 and COVID-19, and financed with REACT-EU resources from the
European regional development fund.

## References

  * [1] K S. Surana.  Transition finite elements for three-dimensional stress analysis.  International Journal for Numerical Methods in Engineering, 15:991–1020, 1980.
  * [2] I. Romero.  Coupling nonlinear beams and continua: Variational principles and finite element approximations.  International Journal for Numerical Methods in Engineering, 114:1192–1212, 2018.
  * [3] I. Steinbrecher, A. Popp, and C. Meier.  Consistent coupling of positions and rotations for embedding 1D cosserat beams into 3d solid volumes.  Computational Mechanics, 69(3):701–732, 2021.
  * [4] A. Furman.  Modeling coupled surface–subsurface flow processes: A review.  Vadose Zone Journal, 7(2):741–756, 2008.
  * [5] M. E. Moghadam, I. E. Vignon-Clementel, R. Figliola, and A. L. Marsden.  A modular numerical method for implicit 0D/3D coupling in cardiovascular finite element simulations.  Journal of Computational Physics, 244:63–79, 2013.
  * [6] C. D’Angelo and A. Quarteroni.  On the coupling of 1d and 3d diffusion-reaction equations: Application to tissue perfusion problems.  Mathematical Models & Methods in Applied Sciences, 18(8):1481–1504, 2008.
  * [7] A. Quarteroni, A. Veneziani, and C. Vergara.  Geometric multiscale modeling of the cardiovascular system, between theory and practice.  Computer Methods in Applied Mechanics and Engineering, 302:193–252, Apr 2016.
  * [8] D. Bernoulli.  Essay d’une nouvelle analyse de la mortalité causée par la petite vérole et des avantages de l’inoculation pour la prévenir.  Mémoires de Mathématiques et de Physique, Académie Royale des Sciences, Paris, pages 1–45, 1760.
  * [9] M. Martcheva.  An Introduction to Mathematical Epidemiology.  Springer Science & Business Media New York, 2015.
  * [10] Malú Grave, Alex Viguerie, Gabriel F. Barros, Alessandro Reali, and Alvaro L. G. A. Coutinho.  Assessing the spatio-temporal spread of covid-19 via compartmental models with diffusion in italy, usa, and brazil.  Archives of Computational Methods in Engineering, 28(6):4205–4223, 2021.
  * [11] Malú Grave, Alex Viguerie, Gabriel F. Barros, Alessandro Reali, Roberto F.S. Andrade, and Alvaro L.G.A. Coutinho.  Modeling nonlocal behavior in epidemics via a reaction–diffusion system incorporating population movement along a network.  Computer Methods in Applied Mechanics and Engineering, 401:115541, 2022.  A Special Issue on Comp. Mod. and Simulation of Infectious Diseases.
  * [12] Alex Viguerie, Alessandro Veneziani, Guillermo Lorenzo, Davide Baroli, Nicole Aretz-Nellesen, Alessia Patton, Thomas E. Yankeelov, Alessandro Reali, Thomas J. R. Hughes, and Ferdinando Auricchio.  Diffusion–reaction compartmental models formulated in a continuum mechanics framework: Application to covid-19, mathematical analysis, and numerical study.  Comput. Mech., 66(5):1131–1152, nov 2020.
  * [13] Nicola Guglielmi, Elisa Iacomini, and Alex Viguerie.  Delay differential equations for the spatially resolved simulation of epidemics with specific application to covid-19.  Mathematical Methods in the Applied Sciences, 45(8):4752–4771, 2022\.
  * [14] F. Laurino and P. Zunino.  Derivation and analysis of coupled PDEs on manifolds with high dimensionality gap arising from topological model reduction.  ESAIM: Mathematical Modelling and Numerical Analysis, 53(6):2047–2080, 2019.
  * [15] M. Kuchta, F. Laurino, K.-A. Mardal, and P. Zunino.  Analysis and approximation of mixed-dimensional pdes on 3D-1D domains coupled with lagrange multipliers.  SIAM Journal on Numerical Analysis, 59(1):558–582, 2021.
  * [16] N. Hagmeyer, M. Mayr, I. Steinbrecher, and A. Popp.  One-way coupled fluid–beam interaction: capturing the effect of embedded slender bodies on global fluid flow and vice versa.  Advanced Modeling and Simulation in Engineering Sciences, 9(1), 2022\.
  * [17] M. Newman.  Networks: An Introduction.  Oxford University Press, 2010.
  * [18] E. Kuhl.  Computational Epidemiology-Data-Driven Modeling of COVID-19.  Springer International Publishing, 2021.
  * [19] J. E. Roberts and J. M. Thomas.  Mixed and hybrid methods.  In P. G. Ciarlet and J. L. Lions, editors, Handbook of Numerical Analysis, pages 523–639. North-Holland, Amsterdam, 1989.
  * [20] F. Brezzi and M. Fortin.  Mixed and hybrid finite element methods.  Springer, Berlin, 1991.
  * [21] F. Auricchio, F. Brezzi, and C. Lovadina.  Mixed finite element methods.  In E. Stein, R. de Borst, and T. J. R. Hughes, editors, Encyclopedia of Computational Mechanics. Wiley & Sons, 2005.
  * [22] R. T. Rockafellar.  Convex analysis, volume 28.  Princeton University Press, Princeton, N.J., 1970.
  * [23] L. C. Evans.  Partial differential equations.  AMS Press, 1999.