# Partial differential equations on hypergraphs and networks of surfaces: derivation and hybrid discretizations

We introduce a general, analytical framework to express and to approximate partial differential equations (PDEs) numerically on graphs and networks of surfaces—generalized by the term hypergraphs. To this end, we consider PDEs on hypergraphs as singular limits of PDEs in networks of thin domains (such as fault planes, pipes, etc.), and we observe that (mixed) hybrid formulations offer useful tools to formulate such PDEs. Thus, our numerical framework is based on hybrid finite element methods (in particular, the class of hybrid discontinuous Galerkin methods).

• 10 publications
• 1 publication
• 9 publications
07/21/2022

### Unsupervised Legendre-Galerkin Neural Network for Stiff Partial Differential Equations

Machine learning methods have been lately used to solve differential equ...
09/28/2020

### Learning Interpretable and Thermodynamically Stable Partial Differential Equations

In this work, we develop a method for learning interpretable and thermod...
08/19/2020

### A Finite Volume Method for Continuum Limit Equations of Nonlocally Interacting Active Chiral Particles

The continuum description of active particle systems is an efficient ins...
11/04/2019

### Abstractions and automated algorithms for mixed domain finite element methods

Mixed dimensional partial differential equations (PDEs) are equations co...
05/14/2020

### Fireshape: a shape optimization toolbox for Firedrake

We introduce Fireshape, an open-source and automated shape optimization ...
06/15/2021

### Structure-preserving Nonlinear Filtering for Continuous and Discontinuous Galerkin Spectral/hp Element Methods

Finite element simulations have been used to solve various partial diffe...
10/26/2020

### An exploratory study on machine learning to couple numerical solutions of partial differential equations

As further progress in the accurate and efficient computation of coupled...

## 1. Introduction

This manuscript establishes a general approach to formulate partial differential equations (PDEs) on networks of (hyper)surfaces, referred to as hypergraphs. Such PDEs consist of differential expressions with respect to all hyperedges (surfaces) and compatibility conditions on the hypernodes (joints, intersections of surfaces). These compatibility conditions ensure conservation properties (in case of continuity equations) or incorporate other properties—motivated by physical or mathematical modeling. We illuminate how to discretize such equations numerically using hybrid discontinuous Galerkin (HDG) methods, which appear to be a natural choice, since they consist of local solvers (encoding the differential expressions on hyperedges) and a global compatibility condition (related to our hypernode conditions). We complement the physically motivated compatibility conditions by a derivation through a singular limit analysis of thinning structures yielding the same results.

Albeit many physical, sociological, engineering, and economic processes have been described by partial differential equations posed on domains which cannot be described as subsets of linear space or smooth manifolds, there is still a lack of mathematical tools and general purpose software specifically addressing the challenges arising from the discretization of these models.

Fractured porous media (see [BerreDK2019] for a comprehensive review) have gained substantial attention and have become an active field of research due to their critical role with respect to flow patterns in several applications in the subsurface, in material science, and in biology. Most commonly, a fracture is described as a very thin, not necessarily planar object in which, for example Darcy’s equation holds. This motivates the singular limit approximation in which a fracture is assumed to be a two dimensional surface within the three dimensional space. When several of these fractures meet, they form a fracture network of two-dimensional surfaces. Thus, fracture networks illustrate a physical application of the type of problem we investigate in this publication. Moreover, a model in which the joints of two (or more) fractures are assigned additional physical properties can be found in [ReichenbergerJBH2006]. Beyond this, fracture networks have been simulated using hybrid high order (HHO) methods [HedinPE2019].

Graph based models for porous media (without fractures) consist of simulating preferential flow paths within the porous matrix. One of the first publications implementing this idea is [Fatt1956] who observed that a network of tubes might approximate the flow of porous media better than the classical model of tube bundles, which has also been used in the most common upscaling techniques—see [SchulzRZRK2019, RayRSK2018] and the references therein for a discussion of those tube models in upscaling procedures. The tube network approach [Fatt1956] has been successfully applied to couple porous media flow to free (Navier–) Stokes flow [WeishauptJH2019].

PDEs on hypergraphs are especially suitable to be used in the description of elastic networks [Eremeyev2019]: Here, we discriminate between one dimensional elastic beam (rod) networks, trusses, etc. and two dimensional elastic plate (shell) networks [LagneseL1993], respectively. Beam networks have been used to model truss bridges and towers (most prominently the Eiffel tower) and other mechanical structures, originating the field of elastic beam theory (for instance [BauchauC2009] for an introduction which also covers elastic plate models). Elastic plate models describe the stability of houses and have several engineering applications such as the description of the stability of (bend) plates (used in automobile industries and several others). They have even been used to understand interseismic surface deformation at subduction zones [KandaS2010].

Elastic beam networks have been used to evaluate elastic constants in amorphous materials. That is, the elastic properties of stiff, beam like polymers have been investigated. Such polymers are key to understanding the cytoskeleton which is an important part of biological cells [Heussinger2007, Lieleg2007], but they are also important for the healing of wounds (fibrin), for skin stability (collagen), and for the properties of paper. Moreover, such models can be used for modelling rubber [PhysRevE.76.031906], foams, and fiber networks [PhysRevLett.96.017802].

Conservation laws in the form of PDEs on hypergraphs have been used in the simulation and optimization of gas networks [RuefflerMH2018] and other networks of pipelines. They have been extended to networks of traffic (streets and data), (tele-)communication, and blood flow. For an overview of the main ideas that are related to these applications, the reader may consult [Garavello2010, BressanCGHP2014]. Additionally, rigorous mathematical analysis of such problems is developing to a field of current research [SkreMR2021].

We conclude the overview over some applications by stating that regular surfaces and volumes can also be interpreted as hypergraphs. Thus, PDEs on surfaces [DziukE2013] and standard “volume” problems (in which the hyperedges have the same dimension as the surrounding space and at most two hyperedges meet in a common hypernode) are also covered by our approach.

Hypergraph models usually are approximations of problems in higher dimensional networks of thin structures, for example a network of thin pipes or thin plates in 3D. As a model example we give a rigorous derivation of a diffusion equation on a hypergraph. More precisely, we consider a network of thin plates in three dimensions, where the thickness of the plates is small compared to their length. We denote the ratio between the thickness and the length by the small parameter

. Due to the different scales the computational effort for numerical simulations is very high. To overcome this problem the idea is to replace the thin-structure by a hypergraph. For this we give a rigorous mathematical justification using asymptotic analysis. We pass to the limit

in the weak formulation of the problem, and derive a limit problem stated on the hypergraph. The solution of this limit-problem is an approximation of the model in the higher-dimensional thin domain. Singular limits for thin plates and shells (leading to lower-dimensional manifolds in the limit ) in elesticity can be found in [ciarlet1997mathematical, ciarlet2000theory]. Dimension reduction for a folded elastic plate is treated in [le1989folded]. Singular limits leading to hypergraphs for fluid equations can be found in [maruvsic2003rigorous], where a Kirchhoff law in a junction of thin pipes is derived, and [maruvsic2019mathematical] where junctions of thin pipes and plates are treated using the method of two-scale convergence.

The remainder of this manuscript is structured as follows: First, we discuss conservation equations on hypergraphs. Second, we rigorously formulate an elliptic model equation and investigate some of its properties in Section 3. Third, we discuss its discretization by means of the HDG method in Section 4. Fourth, we discuss how PDEs on hypergpahs can be obtained by a model reduction approach, in particular, by considering singular limits. The publication is wrapped up, by a section on possible conclusions.

## 2. Conservation equations on geometric hypergraphs

### 2.1. Hypergraphs

A hypergraph consists of a finite set of hyperedges and a finite set of hypernodes. We refer to it as a geometric hypergraph if the hyperedges are smooth, open manifolds of dimension with piecewise smooth, Lipschitz boundary and the hypernodes can be identified with smooth subsets of the boundaries of these hyperedges. More specifically, the boundary of each hyperedge is subdivided into nonoverlapping subsets such that . We associate to

an index vector

and isometries

 ιei:Γei→Nηei,i=1,…,ke. (2.1)

The hypernodes are thus identified with the closures of the subsets of the boundaries of one or more hyperedges. Their dimension is .

has the structure of a hypergraph in the classical sense as each edge connects a set of nodes . The dual hypergraph describes the situation where each hypernode connects hyperedges with indices .

We call a hypernode a boundary hypernode if , i.e., it is part of the boundary of only a single hyperedge. Accordingly, we define the set of boundary hypernodes and the set of interior hypernodes .

As special cases: a geometric graph is a geometric hypergraph where the edges are smooth curves and the nodes are their end points. If every hypernode is either at the boundary or connects exactly two hyperedges, the hypergraph represents a piecewise smooth manifold.

The structure might become more evident if we consider an embedded geometric hypergraph in some ambient space , as in Figure 1 on the left. In this case, the isometries are identical mappings and the hypernodes are identified with the boundary pieces . On the right of this figure, the same hypergraph is displayed without embedding. In this case, the hyperedges are objects in , possibly with a non-flat metric. Hypernodes are intervals in , inheriting their metric through the isometries .

Due to the isometries , every point of a hypernode is uniquely identified with a point on the boundary of each of the hyperedges it connects. Thus, convergence of a point sequence in the union of these hyperedges to a point on the hypernode is well-defined, for instance by considering the (finitely many) subsequences on each hyperedge. Also, a distance between two points on different hyperedges sharing a hypernode is defined locally by these isometries and triangle inequality.

The domain of the hypergraph, its closure, and its boundary are

 ¯¯¯¯Ω=⋃E∈EE∪⋃N∈NN,∂Ω=⋃N∈NBN,Ω=¯¯¯¯Ω∖∂Ω. (2.2)

In this definition, the hyperedges are considered open with respect to their topology and do not contain their boundaries. The hypernodes are closed. We introduce the skeletal domain

 Σ=⋃N∈NN. (2.3)

We make the assumption that is connected. Note that this implies that any two hyperedges are either connected by a common node or not connected, since is open, see (2.2). Without such an assumption, the problems of partial differential equations below separate into subproblems, which then can be analyzed and solved independently.

In Figure 1, comprises all blue hypernodes, which also include the end points of the red hypernode. The union of the red and blue hypernodes is . The domain consists of the interior of the red hypernode and the the three hyperedges.

Many concepts of standard domains in transfer to , even if it is not a manifold. In particular, the notion of a small open ball with radius around , see Figure 2, in is maintained by construction and thus the notion of open subsets. A subset is called compactly embedded in if its closure is contained in and thus has a positive distance to .

A function is continuous on , if it is continuous inside each hyperedge and its limits on a hypernode are consistent between all hyperedges connected by this hypernode. Analogously, a function is in if it is in for all and it is in if it is in for all .

###### Remark 2.1 (Comparison to standard nomenclatures).

In this article, we mix concepts from graph theory, partial differential equations, and finite elements. Thus, a clash of names was unavoidable. What is referred to as a hypernode here, is a face —an edge in two dimensions— in finite element literature, while the hyperedges here correspond to mesh cells or elements. In order to reduce ensuing confusion, we consistently use the term “hyperedge”. Another difference to finite element literature is established by the fact that we consider the hypergraph fixed and are not concerned with refinement limits. Finally, we would like to point out that there has been a concept of geometric hypergraphs in the literature; it is nevertheless very limited, such that we coin this term in a new way here, meaning a hypergraph whose elements are geometric shapes themselves.

### 2.2. Continuity equations on hypergraphs

Next, we conduct a heuristic derivation, employing control volumes

in the shape of infinitesimal, open hyperballs. Let be a conserved quantity and be its flux. Then, the conservation property of is usually stated in integral form such that for any such control volume there holds

 ddt∫Vϱdx=−∫∂VJ⋅ndσ. (2.4)

When is a subset of a Lipschitz manifold , the meaning of this statement is clear if is a smooth tangential vector field in and is the outer normal vector to in the tangential plane of . The term denotes the volume element of the manifold, and is the induced surface element.

If the hyperball intersects a hypernode in which several edges meet, meaning can be given to equation (2.4) by the following observation: if are the hyperedges which meet in inside , then for the intersection has a piecewise smooth boundary . We observe that is in the interior of (see Figure 2 for an illustration) and the boundary of is nowhere tangential to . Thus, with the assumption that no mass is created or destroyed in the hypernode , the conservation property (2.4) can be restated as

 ddt∫Vϱdx=ddtℓ∑i=1∫Viϱdx=−ℓ∑i=1∫∂Vi∖NVJ⋅ndσ. (2.5)

Again, the flux and the outer normal vector to are well defined along in the tangential plane of .

As a generalization of (2.5), we allow for sinks and sources living in the hyperedges and living within hypernode : This can be implemented by setting

 ddt∫Vϱdx=−ℓ∑i=1∫∂Vi∖NVJ⋅ndσ+∫Vfdx+∫NVgdσ, (2.6)

where positive and describe sources, while negative and describes sinks.

Before we convert (2.4) into a problem of partial differential equations, we make the simplifying assumption that the hyperedges and hypernodes are planar and that is the standard Lebesgue measure. This way, we avoid delving into the complexities of surface partial differential equations. This simplification is purely for the ease of presentation and we refer the readers to [DziukE2013] and [BENARTZI2007989] for more general surfaces in the elliptic and hyperbolic settings, respectively.

Thus, in the interior of each hyperedge , we can apply Gauss’ divergence theorem in standard form to obtain

 ∫∂VJ⋅ndσ=∫V∇⋅Jdx. (2.7)

If on the other hand overlaps a hypernode which connects hyperedges , we can still apply the divergence theorem in each hyperedge to obtain

 ∫∂VJ⋅ndσ=ℓ∑i=1[∫Vi∇⋅Jdx−∫NVJ⋅ndσ]. (2.8)

This notion also extends to control volumes intersecting with several hypernodes in a natural way. Then, rearranging the sum over boundary integrals yields

 ∫∂VJ⋅ndσ=∑E∈E∫V∩E∇⋅Jdx−∑N∈NI∫V∩N[[J⋅n]]dσ, (2.9)

for any control volume . Here, is the standard divergence of the differentiable vector field and is the summation operator such that on a hypernode with hyperedges there holds

 [[J⋅n]]=ℓ∑i=1J|Ei⋅nEi. (2.10)

A vector field is usually called solenoidal, if the left side of equation (2.9) vanishes for any control volume . The right hand side of this equation generalizes this notion from standard domains to hypergraphs. Therefore, we call a piecewise smooth vector field solenoidal, if

 ∇⋅J(x) =0 for all x∈E and E∈E, (2.11a) [[J⋅n]](x) =0 for all x∈N and N∈NI. (2.11b)

Note that the second condition is an extension of Kirchhoff’s junction rule from points to higher dimensional hypernodes.

To put it in a nutshell, assuming there are no leaks and sources in hypernodes and hyperedges, the conservation condition (2.5) induces the PDE–interface problem to find and such that

 ∂tϱ+∇⋅J =0 in all E∈E, (2.12a) [[J⋅n]] =0 on all N∈NI. (2.12b)

Analogously, the continuity condition (2.6) induces the PDE interface problem to find and such that

 ∂tϱ+∇⋅J =f in all E∈E, (2.13a) [[J⋅n]] =g on all N∈NI. (2.13b)

In (2.12) and (2.13), and might be linked by some phenomenological description, i.e., (depending on the specific application). Both equations are complemented by appropriate initial and boundary conditions. Beyond this, additional continuity constraints might be formulated, such as , , ….

###### Remark 2.2.

The interface problems in (2.12) and (2.13) resemble the hybrid or hybridized formulation of a PDE, which was introduced for instance in context of the mixed elements of Raviart–Thomas and Brezzi–Douglas–Marini in [RaviartT1977a, RaviartT1977b, BrezziDM1985].

We are, exemplary, going to discuss the continuity equation (2.12) in the context of diffusion problems in section 3.

## 3. Elliptic model equation

The standard diffusion equation in mixed form defined on a hypergraph is a conservation equation of type (2.12) for the flux of a scalar function . This is for instance known as Fourier’s law of thermal conduction, where is the temperature and is the dimensionless heat conductivity of the material. It is also Fick’s law of diffusion where is a concentration and is the diffusion coefficient.

Like in the previous section, we simplify the presentation by assuming that all hyperedges are flat and thus can be identified with a domain in . In the more general case, the differential operators must be replaced by their differential geometric counterparts as in [DziukE2013].

We focus on the stationary case and set the time derivative in (2.12) to zero. Thus, the discussion of the previous section leads to the following problem: find satisfying

 −∇⋅(κ∇u) =f in all E∈E, (3.1a) u =uD on all N∈ND, (3.1b) u|E1 =u|E2 on all N∈N,N⊂∂E1∩∂E2, (3.1c) −[[κ∇u⋅n]] =g on all N∈N∖ND, (3.1d)

for all , right hand sides and , and a diffusion coefficient . A justification by taking the limit of thin domains can be found in Section 5 below.

We observe that in (3.1) the diffusion equation (3.1a) is complemented by three boundary and interface conditions. First, it is closed by a “Dirichlet” boundary condition (3.1b): We choose a non-empty set of “Dirichlet” hypernodes, on which we impose for a prescribed boundary value . In (3.1c), we employ a continuity constraint. This constraint prohibits jumps in the primary unknown across interior nodes, and therefore, loosely speaking, imitates the standard constraint that of the domain.

On interior nodes , we set out with Kirchhoff’s junction law, but with the option of a concentrated source in (3.1d). This equation also incorporates the Neumann condition , since on a boundary hypernode the sum in the definition (2.10) of the operator reduces to a single hyperedge. Note that (3.1d) for on interior nodes serves as compatibility condition for mimicking .

###### Definition 3.1 (Function spaces on hypergraphs).

For each let be the standard Sobolev space on and be the standard trace operator.

Then, we define

 (3.2)

where and are the traces of from the hyperedges and on , respectively. Due to the equality of traces in the definition of , we can define the trace operator to the skeleton

 γ:H→M:={μ∈L2(Σ)∣∣∣μ|∂E∈H1/2(∂E)for all E∈E}. (3.3)

Additionally, the spaces and are defined as

 M0:= {μ∈M:μ|N=0 for all N∈ND}, (3.4) H0:= {u∈H:γu∈M0}, (3.5)

and we denote the dual spaces of by and of by .

Norms ( and ) on the respective spaces ( and ) are induced by summed versions of the local scalar-products:

 (u,v)H:= ∑E∈E(u|E,v|E)H1(E), ∥u∥2H:= (u,u)H, (3.6) ⟨λ,μ⟩M:= ∑E∈E⟨λ|E,μ|E⟩H1/2(∂E), ∥μ∥2M:= ⟨μ,μ⟩M, (3.7)

such that

 ∥u∥2H=∑E∈E∥u|E∥2H1(E)and∥μ∥2M=∑E∈E∥μ|∂E∥2H1/2(∂E). (3.8)

These definitions have a few immediate consequences:

1. is a well-defined and surjective, linear, and continuous operator.

2. We have the Gelfand triple relations

 H0↪L2(Ω)≅[L2(Ω)]⋆↪H⋆0, (3.9) M0↪L2(Σ)≅[L2(Σ)]⋆↪M⋆0. (3.10)

Note that is analogous to space in of Raviart and Thomas [RaviartT1977a].

###### Lemma 3.2.

The space with inner product is a Hilbert space.

###### Proof.

Obviously, is a subspace of the Hilbert space , and the function

 h:⨁E∈EH1(E)∋u↦∑N∈NI∑¯E1,¯E2⊃N∥u|E1−u|E2∥2L2(N)∈R

is continuous and is its kernel. Thus, is closed. ∎

###### Definition 3.3.

A weak solution to the primal formulation of (3.1) with , , and is a function with on all , and

 ∑E∈E∫Eκ∇u⋅∇vdx=⟨f,v⟩H⋆0,H0−⟨g,v⟩M⋆0,M0∀v∈H0. (3.11)

In particular, if and , we can rewrite (3.11) as

 ∑E∈E∫Eκ∇u⋅∇vdx=∫Ωfvdx−∑N∈N∫Ngvdσ∀v∈H0. (3.12)

### 3.1. Existence and uniqueness of solutions

###### Theorem 3.4.

Assume for that there is a lifting with on all . If with a. e., , , and all are Lipschitz domains, there is an unique weak solution according to Definition 3.3, which continuously depends on the data.

###### Proof.

Due to the existence of , we can reduce the problem to the one with homogeneous Dirichlet values if we replace by and modifying the right hand side accordingly. Since the right hand side is bounded and is a Hilbert space, it suffices to show ellipticity of the weak form to conclude the proof by the Lax–Milgram lemma. We note that for there holds

 ∑E∈E∫Eκ∇v⋅∇vdx≥κ0∑E∈E∥∇v∥2L2(E) (3.13)

Thus, the following Poincaré–Friedrichs inequality implies ellipticity and concludes the proof. ∎

###### Lemma 3.5 (Poincaré–Friedrichs inequality for H0).

For all it holds that

 ∥v∥L2(Ω)≤C∑E∈E∥∇v∥L2(E).
###### Proof.

Similar to the standard case of subdomains in , this inequality follows easily by contradiction: To this end, we assume that there is a sequence with

 ∥vn∥L2(Ω)=1and∑E∈E∥∇vn∥L2(E)≤1n. (3.14)

Thus, is bounded in for all . Hence, by the weak compactness of the unit ball in and the Rellich-Kondrachov theorem, there exists a subsequence (also denoted ) such that

 vn|E→~v|E in L2(E),andvn|E⇀~v|E in H1(E). (3.15)

We have that (due to its completeness), and that the seminorm . Thus, is constant in all and overall continuous. Therefore it is overall constant and has to be zero, due to the zero boundary condition on Dirichlet nodes and the connectedness of . Hence, the strong convergence of in implies , which contradicts . Therefore, the Poincaré–Friedrichs inequality is valid. ∎

## 4. HDG method for elliptic model equation

When we derived PDE problems on hypergraphs, we were led to a formulation local on each hyperedge with coupling conditions on hypernodes. This is a structure which is nicely reflected in hybridized methods. Indeed, there the separation goes one step further. By putting degrees of freedom on the hypernode, values on hyperedges are not coupling anymore to other hyperedges across these hypernodes, but only to the values on the hypernodes constituting their boundary. Thus, differing from standard or discontinuous finite element methods, the number of hyperedges attached to a hypernode does not affect the solution process on a single hyperedge. Therefore, we consider hybridized methods ideally suited to PDEs on hypergraphs.

Hybridized discontinuous Galerkin (HDG) methods break the continuity condition (3.1c) by introducing Lagrange multipliers on each hypernode which enforce the continuity of fluxes (3.1d) weakly. It turns out though, that the Lagrange multiplier is an approximation to the solution of (3.1) on the skeleton itself.

With such methods, the actual PDE (3.1a) is represented locally on each hyperedge by Steklov-Poincaré operators on the hyperedges, which transform function values to flux values on the boundary of the hyperedges, a process called “local solver” in HDG terminology. The global problem is posed in terms of the degrees of freedom on the hypernodes only, yielding a square, linear system of equations.

In this respect, HDG methods have a similar structure as the family of HHO methods. These are based on hybridizing the primal formulation and lead to a rather simple error analysis on polytopic meshes where only projections are used (as opposed to the rather complicated projections used for HDG). This is achieved by a novel stabilization design [DiPietro2015]. For recent developments in hybrid high-order and HDG methods, the reader may consult [Burman2018, Qiu2016].

The separation of the local solution of bulk problems from the global coupling of interface variables is also achieved by the virtual element method [VEM0, VEM1]. Thus, it fits into our view of coupled differential equations on connected hyperedges. Different to the methods discussed so far, it does not rely on polynomial shape functions inside mesh cells but rather on forms of fundamental solutions of any shape [VEM2]. Accordingly, when applied to hypergraphs, the actual type of local solvers and of the specific boundary trace operators will differ from our approach, but remain within the same principal concept.

### 4.1. The hybridized dual mixed formulation

In physical applications, there often is a need to receive reasonable approximations for both the primal unknown and the dual unknown . In other words, considering diffusion, we would like to know both the distribution of some species’ concentration and the species “movement” (flux). This becomes particularly important if we interpret as pressure and as fluid flow through a porous medium (Darcy’s equation). In this situation, the flow field will govern the movement of chemical species dissolved within the fluid. It is often the main quantity of interest and conservativity is crucial. Therefore, we turn to the mixed formulation.

The mixed HDG methods use the weak, dual, mixed, hybrid formulation of (3.1), i.e., find with on all such that

 ∫E[u(∇⋅p)−κ−1q⋅p]dx =∫∂Eλp⋅ndσ ∀p∈⨁E∈EHdiv(E), (4.1a) ∫Ev∇⋅qdx =∫Efvdx ∀v∈L2(Ω), (4.1b) ∑E∈E∫∂E(q⋅n)μdσ =⟨g,μ⟩M⋆0,M0 ∀μ∈M0. (4.1c)

Well-posedness of this formulation can be deduced from Theorem 3.4 if . Indeed, on the one hand, this implies that the (uniquely existing) solution of Definition 3.3 solves (4.1) with , , and . On the other hand, for any solution of (4.1), we have (by the space’s definition), and in the weak sense. Therefore, any solution to (4.1) satisfies Definition 3.3.

Equations (4.1a) & (4.1b) are local equations on the hyperedge, like in the standard case of a domain. They only couple to the Lagrange multipliers on the boundary of the hyperedge. Thus, we can eliminate them locally in the fashion of the Schur complement method. To this effect, we introduce the local solution operator for the right hand side . It is in fact a Steklov-Poincaré operator on mapping the Dirichlet data to the normal trace of the flux in (4.1c).

Then, the solution of (4.1) can be characterized as

 ∑E⟨SEλ,μ⟩M⋆0,M0=⟨g,μ⟩M⋆0,M0. (4.2)
###### Remark 4.1.

The Steklov-Poincaré operators in this equation are the same ones as in the case of a manifold. They do not depend on the connectivity of a hypernode to other hyperedges. Thus, their implementation does not differ from that of a standard finite element method. The only difference lies in the structure of the sum on the left, and is thus almost purely of algebraic nature.

For inhomogeneous right hand side , we can define the operators

 (~U,~Q):L2(Ω)→L2(Ω)×⨁E∈EHdiv(E)f↦(u,q) (4.3)

by an edge-wise solution of (4.1a) and (4.1b) with .

In order to obtain a better understanding of the Steklov-Poincaré operators, we follow the route laid out in [CockburnGL2009] for the discrete version and define the solution operators

 (4.4)

which map a given to the element-wise solution of (4.1a) and (4.1b) with .

The well-posedness and linearity of all local solution operators follow directly from the fact (see for instance [BoffiBrezziFortin13]) that the mixed formulation on a single hyperedge is well-posed for any given and its solution depends continuously on . Entering these solution operators into (4.1c) yields

 −∑E∈E∫∂E(Qλ⋅n)μdσ=∑E∈E∫∂E(~Qf⋅n)μdσ−⟨g,μ⟩M⋆0,M0, (4.5)

since . By some simple transformations of (4.1), we can write (4.2) with inhomogeneous right hand sides in terms of bilinear and linear forms. This argument allows to reduce the problem to finding with on all , such that

 a(λ,μ)=b(μ)∀μ∈M0, (4.6a) a(λ,μ) =∑E∈E∫Eκ−1Qλ⋅Qμdx, (4.6b) b(μ) =∑E∈E∫∂E(~Qf⋅n)μdσ−⟨g,μ⟩M⋆0,M0. (4.6c)

Obviously, bilinear form and linear form are continuous due to the continuity of operators and . Surprisingly, we have recovered a symmetric bilinear form. The following lemma is a key to the discrete well-posedness and adds the fact that this form is even -elliptic.

###### Lemma 4.2.

If and , bilinear form from (4.6b) is elliptic.

###### Proof.

Like in the proof of the Poincaré-Friedrichs inequality, we prove ellipticity of by a contradiction argument. To this end, let a sequence in such that

 ∥λn∥M=1and∥Qλn∥L2(Ω)→0. (4.7)

Thus, there exists a subsequence in , and by the compact embedding of in there holds again for a subsequence in . Since is continuous and converges weakly in we obtain weakly in for each hyperdege and therefore also in . Thus, (4.7) implies that . We denote by the solution to (4.1a) and (4.1b) associated to , especially we have . Hence, for every we have

 ∫Eu(∇⋅p)dx=∫∂Eλp⋅ndσ∀p∈Hdiv(E). (4.8)

Moreover, we have that the divergence operator

 div:H10(E)d→L20(E):={u∈L2(E):∫Eu=0} (4.9)

is surjective, and therefore

 ∫Euφdx=0∀φ∈L20(E). (4.10)

This, however, implies that is constant on . Furthermore, since is constant on , we can deduce by (4.8) and Gauss’ divergence theorem that is constant and

 ~λ≡conston∂E∀E∈E. (4.11)

The contradiction argument is concluded by the fact that some hyperedges are adjacent to the Dirichlet nodes and thus on their boundary. For the other hyperedges, follows from connectedness of . Thus, on in contradiction to . Altogether we showed that for all it holds that

 ∥λ∥M≲∥Qλ∥L2(Ω).

Together with (4.6b) we obtain for a positive constant

 a(λ,λ)≥α∥λ∥2M,

i.e., the ellipticity of on .

Thus, the Lax-Milgram lemma guarantees a unique solution as soon as it is clear that and generate a right hand side in the dual of . But for this is obvious as . For , there is a unique solution of (4.1a) and (4.1b) with . Its trace is in for each hyperedge and thus bounded on . ∎

### 4.2. HDG methods in dual mixed form

Let be some finite dimensional, scalar function space. Then, we define the space of discrete functions on the skeleton by

 M:={λ∈L2(Σ)∣∣ ∣∣λ|N∈^M∀N∈Nλ|N=0∀N∈ND}. (4.12)

The mixed HDG methods involve a local solver on each hyperedge , producing hyperedge-wise approximations and and of the functions and in equation (4.1), respectively. Here, is some finite dimensional, scalar function space, and is some finite dimensional, vector valued function space. We will also use the concatenations of the spaces and , respectively, as a function space on , namely

 V:={v∈L2(Ω)∣∣v|E∈VE,∀E∈E},W:={q∈L2(Ω;Rd)∣∣q|E∈WE,∀E∈E}. (4.13)

The HDG scheme for (4.1) on a hypergraph consists of the local solver and a global coupling equation. The local solver is defined hyperedge-wise by a weak formulation of (4.1) in the discrete spaces and defining suitable numerical traces and fluxes. Namely, given find and , such that

 ∫E1κQE⋅pdx−∫EUE∇E⋅pdx =−∫∂Eλp⋅ndσ (4.14a) ∫∂E(QE⋅n+τUE)vdσ−∫EQE⋅∇Evdx =τ∫∂Eλvdσ (4.14b)

hold for all , and all , and for all . Here, is the penalty coefficient. While the local solvers are implemented hyperedge by hyperedge, it is helpful for the analysis to combine them by concatenation. Thus, the local solvers define a mapping

 M→V×Wλ↦(Uλ,Qλ), (4.15)

where for each hyperedge holds and . Analogously, we set and , where now the local solutions are defined by the system

 ∫E1κQE⋅pdx−∫EUE∇E⋅pdx =−∫∂EuDp⋅ndσ (4.16a) ∫∂E(QE⋅n+τUE)vdσ−∫EQE⋅∇Ev dx (4.16b) = ∫Efvdx+τ∫∂EuDvdσ

Once has been computed, the HDG approximation to (4.1) on will be computed as

 U=Uλ+U(f,uD),Q=Qλ+Q(f,uD) (4.17)

The global coupling condition is derived through a discontinuous Galerkin version of mass balance and reads: Find , such that for all

 (4.18)
###### Remark 4.3.

Hybridized DG methods in dual mixed form differ by the choice of local polynomial spaces and the stabilization parameter . Defining as the space of multivariate polynomials of degree at most , Table 1 lists some well-known combinations on simplices.

Well-posedness of the local solvers for all of them is proven in [CockburnGL2009]

and the works cited there. Analogous combinations based on tensor product polynomials exist for hypercubes.

Existence and uniqueness of the discrete solution , , and to the HDG method can be shown repeating the arguments mentioned in Section 4.1 in the finite-dimensional setting. A natural assumption is the well-posedness of the local problems (4.14), see Remark 4.3.

Given the local solvers, the HDG method for elliptic diffusion problems is consistent with respect to the solution to (4.1). Using consistency, we can immediately apply the analysis in [CockburnGW2009], as it proceeds locally for each hyperedge. Thus, we obtain optimal convergence rates for LDG-H (and also RT-H by slight adaptions) on simplicial hypergraphs. They also transfer to quadrilateral hypergraphs, since these allow for a Raviart–Thomas projection satisfying equation (2.7) in [CockburnGW2009].

### 4.3. Numerical convergence tests for LDG-H

Next, we consider a convergence example on a hypergraph. It is constructed to approximate

 −∇⋅(κ∇u)=f in E∈E,u=uD%  on N∈ND, (4.19)

where the Dirichlet nodes are those that are located on the boundary of with .

The filling indicates that the cube has been times uniformly refined (in the standard three dimensional sense), and the calculation is conducted on the dimensional “surfaces” of this filling. These surfaces themselves might be further refined times, and these refined surfaces are identified to be our standard hyperedges, see Figure 3 for an illustration.

The dimensional faces of this approach are interpreted as nodes and the nodes located on the boundary of the unit cube are considered Dirichlet nodes. All other nodes are supposed to be interior nodes. The solution is constructed to be , diffusion coefficient , and right-hand side . Of course, polynomial degrees are supposed to exactly reproduce the given solution, which is true in our numerical experiments. Thus, we only plot the errors for in Table 2.

Interestingly, the errors converge although with filling , also the computational domain increases for and . However, the rate of convergence deteriorates by if , and if . The optimal order is obtained for .

Beyond this, the refinement indicated with uses filling level and then uniformly refines the respective faces. This does not lead to an increase of the computational domain (even if ) and, therefore, gives the optimal convergence rate .

The aforementioned results have been obtained using our code HyperHDG [HyperHDGgithub].

## 5. Hypergraph PDE as singular limit

The aim of this section is to derive the hypergraph model (3.1) as a singular limit of a 3D-model problem as illustrated in Figure 4. We exemplary use the figure to illuminate the basic ideas: We assume to have a diffusion problem on a domain consisting of three thin plates (in gray) and a (red) joint. This is the problem, which we would like to solve. However, we do not want to solve it in three spatial dimensions, but would like to reduce it to a two-dimensional problem—for example since we have limited compute sources, the domain is very large, or the domain is very complicated to mesh. Thus, we let the thickness of the three plates (and therefore also the thickness of the joint) go to zero by considering , and construct a two dimensional limit problem. The solution of this two dimensional problem lives on the mid-planes of the three planes and their joint. It in some sense is supposed to approximate the solution of the original (three-dimensional) problem for which is a small, positive number.

The principal idea of the limit process is to map equations on the thin structures depending on to fixed reference domains independent of , where we can use standard compactness methods from functional analysis. However, the transformed problem includes -dependent coefficients. Thus, the crucial point for the derivation of the limit model is to establish a priori estimates that are uniform with respect to .

### 5.1. Description of the 3D model problem

We consider the simplified case of one hypernode with length connecting hyperedges for which are rectangles with side lengths and . Thus, Figure 4 shows the case with . The opposite node of with respect to is denoted by (and is a boundary node). Without loss of generality, we assume that lies in the -axis and we have

 N={se1:s∈(0,L)}. (5.1)

We denote by a unit normal vector to and define extruded hyperedges for and

 ~Eϵi:={x∈R3:x=y+sνi for% s∈(−diϵ2,diϵ2),y∈Ei}.

Hence, is a hexahedron with side lengths and , and with thickness . We construct now a domain which contains the union of all these extruded hyperedges and a nonoverlapping decomposition of this domain. To this end, let be chosen such that the sets

 Eϵi:={x∈~Eϵi:dist(x,~Sϵi)>αϵ}, (5.2)

do not overlap. We denote the side of that contains by , and define

 Sϵi :=int{x∈∂Eϵi:dist(x,~Sϵi)=αϵ}.

The side lengths of are and . Additionally, we define the convex hull of the node and the sides :

 Nϵ:=int(conv{¯¯¯¯¯N,¯¯¯¯¯¯Sϵ1,…,¯¯¯¯¯¯¯Sϵm}).

By construction, we have

 Nϵ=(0,L)×ϵω.

Then, we define the thin domain as

 Ωϵ:=Nϵ∪m⋃i=1(Eϵi∪Sϵi).

On we define a diffusion problem and then to pass to the limit in order to derive a problem on the hypergraph . To ensure uniqueness for our model we assume a zero-Dirichlet boundary condition on one face . We consider the following problem for the unknown

 −∇⋅(κϵ∇uϵ) =fϵ in Ωϵ, (5.3a) uϵ =0 on SϵD, (5.3b) −κϵ∇uϵ⋅ν =