# Use of algebraic dual representations in domain decomposition methods for
Darcy flow in 3D domains

V. Jain<EMAIL_ADDRESS>A. Palha  M. Gerritsma  Delft University of
Technology, Faculty of Aerospace Engineering, P.O. Box 5058, 2600 GB Delft,
The Netherlands

###### Abstract

In this work we use algebraic dual representations in conjunction with domain
decomposition methods for Darcy equations. We define the broken Sobolev spaces
and their finite dimensional counterparts. In addition, a global trace space
is defined that connects the solution between the broken spaces. Use of dual
representations results in a sparse metric free representation of the
constraint on divergence of velocity, the pressure gradient term and on the
continuity constraint across the sub domains. To demonstrate this, we solve
two test cases: i) manufactured solution case, and ii) industrial benchmark
reservoir modelling problem SPE10. The results demonstrate that the domain
decomposition scheme, although with more unkowns, requires less memory and
simulation time as compared to the continuous Galerkin formulation.

###### keywords:

Domain decomposition methods, Algebraic dual spaces, Darcy equations, SPE10,
Mimetic spectral element method, Hybrid finite elements

##  1 Introduction

Since the early days of computational fluid dynamics there has been an immense
increment in computational power. However, the need for development of fast
numerical algorithms has persisted consistently. To that motive, the class of
domain decomposition (DD) methods has played an important role in reducing the
run times and memory requirements for numerical simulations. The general aim
for these methods is that the domain of the problem is broken up into a set of
many small independent sub domains. New degrees of freedom are introduced, the
Lagrange multipliers, that impose the appropriate continuity constraints
across these domains. It is then possible to solve for the coupled system of
Lagrange multiplier equations only, which is a smaller system compared to the
continuous unbroken formulation. The local solution in the sub domains is
obtained independently of each other as a function of the solution of the
Lagrange multipliers. This process reduces the computational burden of the
problem in the sense that it is less demanding in terms of required memory and
computational time, without compromising on the accuracy of the results. Over
the last decades, many variants of DD methods have been developed, for
example, hybrid method [1], mortar method [2], FETI method [3, 4], DPG method
[5, 6, 7], Steklov-Poincaré method [8, 9], etc. For a comprehensive discussion
of these methods we refer the reader to the review papers [10, 11]. The two
main challenges for DD methods are: i) to find a subset of suitable finite
dimensional spaces such that the system matrices are not singular, i.e. they
do not produce spurious kinematic modes [12], especially with respect to the
trace spaces of the Lagrange multipliers; and ii) how to efficiently solve the
global system of Lagrange multiplier equations which can become large for
practical applications and is characterized by high condition numbers. The
primary objective of this paper is to address the first challenge, i.e. to
define the framework and extend the use of algebraic dual representations
introduced in [13] for DD formulation of Darcy flow. The DD formulation used
in this work is based on the hybrid method which is a form of discontinuous
Galerking formulation. It is shown that the use of algebraic dual
representation results in a sparse metric free representation for the
divergence constraint on the velocity, the pressure gradient and the
continuity constraint across the sub domain, even for high order spectral
element method. We also see that the continuity constraints are local to the
boundary face of adjoining sub domain elements. This construction is similar
to the use of dual spaces in [2], where they have non-confirming sub domain
boundaries and the continuity constraints are also local to the elements on
the adjoining sub domain boundaries only.

In this paper we use spectral elements of order $N=1,2,3$ and demonstrate the
advantage of using DD formulation with algebraic dual representations using
two test cases. The first test case is a manufactured solution taken from [14]
on a randomly deformed smooth domain. Using this case we first show that using
the DD formulation gives the same solution as the continuous formulation and
have speed-up in simulation times. It is observed that on mesh refinement the
speed-up in simulation time is larger in the case of high order elements (see
Table 3). We also show optimal rates of convergence for varying mesh
refinements. The second test case is a benchmark reservoir modelling test case
SPE10 [15] that demonstrates the use of this formulation on a practical
application. The purpose of this work is fourfold: i) to demonstrate that the
algebraic dual representations can be extended to the framework of DD methods
without compromising on the accuracy of results, ii) to show that using DD
method we can go to higher levels of mesh refinement with same memory
hardware, iii) the use of DD formulation is more efficient in terms of
simulation time, and iv) to show the applicability of the method for industry
benchmark problems such as the SPE10 case.

This paper is structured as follows: In Section 2 we define the broken Sobolev
spaces for DD formulation. In Section 3 we define the finite dimensional
subset of the broken Sobolev spaces. In Section 4 we state the Darcy problem
and the weak formulation of this problem for DD method. The algebraic
formulation for this problem and the solution steps are also described in this
section. In Section 5 we present the results for the two test cases, i) the
manufactured test case from [14], and ii) the benchmark test case SPE10 [15].
We draw conclusions and discuss the scope for future work in Section 6.

##  2 Broken Sobolev spaces

Let $\Omega\subset\mathbb{R}^{3}$ be a bounded domain with Lipschitz boundary
$\partial\Omega$. Let $L^{2}\left(\Omega\right)$ be the space of square
integrable functions and $[L^{2}\left(\Omega\right)]^{3}$ the space of square
integrable vector fields in 3D, and $H^{1}\left(\Omega\right)$,
$H\left(\mathrm{div};\Omega\right)$ be the Hilbert spaces defined as

$\begin{array}[]{ll}H^{1}\left(\Omega\right)&:=\left\\{p\in
L^{2}\left(\Omega\right):\mathrm{grad}\
p\in[L^{2}\left(\Omega\right)]^{3}\right\\}\\\\[4.73611pt]
H\left(\mathrm{div};\Omega\right)&:=\left\\{\bm{v}\in[L^{2}\left(\Omega\right)]^{3}:\mathrm{div}\
\bm{v}\in L^{2}\left(\Omega\right)\right\\}\end{array}\;.$

The trace spaces of $H^{1}\left(\Omega\right)$ is defined as

$\begin{array}[]{ll}H^{1/2}\left(\partial\Omega\right)&:=\left\\{\lambda\in
L^{2}\left(\partial\Omega\right):\exists\ p\in
H^{1}\left(\Omega\right)\quad\mathrm{s.t.}\quad\lambda=p|_{\partial\Omega}\right\\}\\\\[4.73611pt]
\end{array}\;,$

and we denote by $H^{-1/2}\left(\partial\Omega\right)$ its dual space.

Let $\Omega$ be broken into $T$ non-overlapping open sub domains
$\mathcal{M}_{i}$ with Lipschitz boundary $\partial\mathcal{M}_{i}$,
$i=1\;,\ldots\;,T$, such that

$\Omega=\bigcup_{i=1,\ldots\;,T}\mathcal{M}_{i}\qquad\mbox{with,}\qquad\mathcal{M}_{i}\subset\Omega\qquad\mbox{and}\qquad\overset{\circ}{\mathcal{M}}\
_{i}\bigcap\overset{\circ}{\mathcal{M}}\ _{j}=\emptyset\quad\mbox{for}\quad
i\neq j\;.$ (1)

Let $\Omega_{T}$ be the set of sub domains, and $\partial\Omega_{T}$ be the
set of boundaries of these sub domains defined as

$\Omega_{T}=\left\\{\mathcal{M}_{i}\right\\}_{i=1,\ldots\;,T}\qquad\qquad\partial\Omega_{T}=\left\\{\partial\mathcal{M}_{i}\right\\}_{i=1,\ldots\;,T}\;.$

We define the broken Sobolev spaces for the set of sub domains $\Omega_{T}$ as

$\displaystyle L^{2}\left(\Omega_{T}\right):=$
$\displaystyle\prod_{\mathcal{M}\in\Omega_{T}}L^{2}\left(\mathcal{M}\right)$
$\displaystyle H\left(\mathrm{div};\Omega_{T}\right):=$
$\displaystyle\prod_{\mathcal{M}\in\Omega_{T}}H\left(\mathrm{div\;;\mathcal{M}}\right)\;.$
$\displaystyle H^{-1/2}\left(\partial\Omega_{T}\right):=$
$\displaystyle\prod_{\Gamma\in\partial\Omega_{T}}H^{-1/2}\left(\Gamma\right)$

Let the set of boundary faces, $\partial\Omega_{e}$ be defined as

$\partial\Omega_{e}=\left\\{\gamma_{ij}=\partial\overline{\mathcal{M}}_{i}\bigcap\partial\overline{\mathcal{M}}_{j}\qquad\mbox{for}\qquad
i<j\right\\}\;.$
$\partial\Omega_{\Gamma}=\left\\{\gamma_{i}=\partial\mathcal{M}_{i}\bigcap\partial\Omega\qquad\mbox{for}\qquad
i=1,\ldots,T\right\\}\;.$

Then we define the trace space

$H^{1/2}\left(\partial\Omega_{e}\right)=\left\\{\lambda|_{\gamma}\in
L^{2}\left(\gamma\right),\gamma\in\partial\Omega_{e}\right\\}=\prod_{\Gamma\in\partial\Omega_{e}}H^{1/2}\left(\Gamma\right)\;.$

### Notation

We will denote the $L^{2}$-inner product by $\left(\cdot\;,\cdot\right)$. A
set of coordinates is denoted by $\bm{x}$. We denote by $N$ the highest
polynomial degree of basis functions used in an element is $N-1$.

We use two set of finite element spaces, the primal representation and the
dual representation from [13]. The primal representation of finite element
space, basis functions, and degrees of freedom are dentoed by
$X\left(\cdot\right)$, $\Psi\left(\cdot\right)$, $\mathcal{N}$ respectively.
The dual representation of space $X\left(\cdot\right)$, its associated basis
functions and expansion coefficients are denoted with a tilde as
$\widetilde{X}\left(\cdot\right)$, $\widetilde{\Psi}\left(\cdot\right)$ and
$\widetilde{\mathcal{N}}\left(\cdot\right)$, respectively. The basis functions
are always represented as row vectors

$\Psi^{k}\left(\bm{x}\right)=\left(\bm{\epsilon}^{(k)}_{1}\left(\bm{x}\right)\quad\bm{\epsilon}^{(k)}_{2}\left(\bm{x}\right)\quad\ldots\quad\bm{\epsilon}^{(k)}_{d_{k}}\left(\bm{x}\right)\right)\;,$

where $d_{k}$ is the dimension of polynomial vector space of
$X\left(\cdot\right)$. The expansion coefficients are always represented as
column vectors

$\mathcal{N}^{k}\left(\bm{u}\right)=\left(\begin{array}[]{c}\mathcal{N}_{1}^{k}\left(\bm{u}\right)\\\\[4.73611pt]
\mathcal{N}_{2}^{k}\left(\bm{u}\right)\\\\[4.73611pt] \vdots\\\\[4.73611pt]
\mathcal{N}_{d_{k}}^{k}\left(\bm{u}\right)\end{array}\right)$

We will Lemma 2 from [13], which states

######  Lemma 1.

If $\Psi^{k}\left(\bm{x}\right)$ and $\widetilde{\Psi}^{k}\left(\bm{x}\right)$
are basis functions from primal and dual representations respectively, then
these bases are bi-orthogonal with respect to each other

$\int_{\Omega}\widetilde{\Psi}^{3-k}\left(\bm{x}\right)\Psi^{k}\left(\bm{x}\right)\mathrm{d}\Omega=\mathbb{I}\;,$
(2)

where $\mathbb{I}$ is the identity matrix of dimension $d_{k}$.

######  Corollary 1.

The inner product between variables of primal and dual representation is
vector dot product of expansion coefficients.

###### Proof.

Let $p,q\in X\left(\Omega\right)\times\widetilde{X}\left(\Omega\right)$, the
the inner product is given by

$\int_{\Omega}p\ q\
\mathrm{d}\Omega=\widetilde{\mathcal{N}}^{3-k}\left(q\right)^{\intercal}\left(\int_{\Omega}\widetilde{\Psi}^{3-k}\left(\bm{x}\right)^{\intercal}\Psi^{k}\left(\bm{x}\right)\mathrm{d}\Omega\right)\mathcal{N}^{k}\left(p\right)\stackrel{{\scriptstyle\eqref{eq:primal_dual_int}}}{{=}}\widetilde{\mathcal{N}}^{3-k}\left(q\right)^{\intercal}\mathcal{N}^{k}\left(p\right)\;.$

∎

In general, if not explicitly mentioned otherwise, we use Gauss-Lobatto-
Legendre points for numerical integration.

##  3 Finite dimensional spaces

In this section we will define the finite dimensional spaces for hexahedral
elements in 3D domains. The domain $\Omega$ is discretized into a mesh that
consists of points, edges, surfaces and volumes. For $N>1$ each element also
consists of a GLL mesh. We will first introduce the finite dimensional spaces
for the domain $\Omega$ as defined in [13]. The basis functions for domain
$\Omega$ are obtained from mapping of basis functions on a reference domain
$\widehat{\Omega}=[-1,1]^{3}$. For further details on construction we refer
the reader to [13, §4.5]. Here we directly use the mapped basis functions and
the relevant properties of these spaces. We will use the primal
representations to define finite element spaces: i)
$D\left(\Omega\right)\subset H\left(\mathrm{div};\Omega\right)$, ii)
$S\left(\Omega\right)\subset L^{2}\left(\Omega\right)$, iii)
$D_{b}\left(\partial\Omega\right)\subset H^{-1/2}\left(\partial\Omega\right)$,
and dual representations to define the finite element spaces: i)
$\widetilde{S}\left(\Omega\right)\subset L^{2}\left(\Omega\right)$, ii)
$\widetilde{D}_{b}\left(\partial\Omega\right)\subset
H^{1/2}\left(\partial\Omega\right)$.

These spaces will then be used to define the broken finite dimensional spaces
that will be used for DD formulation.

###  3.1 Primal representations

####  3.1.1 Finite element space $D\left(\Omega\right)\subset
H\left(\mathrm{div};\Omega\right)$

Let $\bm{u}$ be the flux component of a vector field. For an element
$\bm{u}\in D\left(\Omega\right)$, let $\Psi^{2}\left(\bm{x}\right)$,
$\mathcal{N}^{2}\left(\bm{u}\right)$ form the row vector of basis functions
and the column vector of expansion coefficients, given as

$\Psi^{2}\left(\bm{x}\right)=\left(\bm{\epsilon}^{(2)}_{1}\left(\bm{x}\right)\quad\bm{\epsilon}^{(2)}_{2}\left(\bm{x}\right)\quad\ldots\quad\bm{\epsilon}^{(2)}_{d_{2}}\left(\bm{x}\right)\right)\;,\qquad\qquad\mathcal{N}^{2}\left(\bm{u}\right)=\left(\begin{array}[]{c}\mathcal{N}_{1}^{2}\left(\bm{u}\right)\\\\[4.73611pt]
\mathcal{N}_{2}^{2}\left(\bm{u}\right)\\\\[4.73611pt] \vdots\\\\[4.73611pt]
\mathcal{N}_{d_{2}}^{2}\left(\bm{u}\right)\end{array}\right)\;,\qquad\qquad\mbox{respectively.}$
(3)

Here, $\mathcal{N}^{2}_{i}\left(\bm{u}\right)=\int_{f_{i}}\bm{u}\cdot\bm{n}\
\mathrm{dA}$ denotes the net flux $\bm{u}$ through the face $f_{i}$.

Then we can represent $\bm{u}$ as

$\bm{u}\left(\bm{x}\right)=\Psi^{2}\left(\bm{x}\right)\mathcal{N}^{2}\left(\bm{u}\right)\;.$

Using this notation, for two elements $\bm{u},\bm{v}\in D\left(\Omega\right)$
the $L^{2}$-inner product is given as

$\left(\bm{u}\;,\bm{v}\right)=\int_{\Omega}\bm{u}^{\intercal}\bm{v}\mathrm{d}\Omega={\mathcal{N}^{2}\left(\bm{u}\right)}^{\intercal}\left(\int_{\Omega}{\Psi^{2}\left(\bm{x}\right)}^{\intercal}\Psi^{2}\left(\bm{x}\right)\mathrm{d}\Omega\right)\
\mathcal{N}^{2}\left(\bm{v}\right)={\mathcal{N}^{2}\left(\bm{u}\right)}^{\intercal}\mathbb{M}^{(2)}\mathcal{N}^{2}\left(\bm{v}\right)\;,$

where,
$\mathbb{M}^{(2)}:=\int_{\Omega}{\Psi^{2}\left(\bm{x}\right)}^{\intercal}\Psi^{2}\left(\bm{x}\right)\mathrm{d}\Omega$
is the mass matrix associated with the basis functions
$\Psi^{2}\left(\bm{x}\right)$.

If there is a symmetric positive definite permeability tensor $\mathbb{K}$,
then the weighted inner product is given as

$\left(\bm{u}\;,\bm{v}\right)_{\mathbb{K}^{-1}}=\int_{\Omega}\bm{u}^{\intercal}\mathbb{K}^{-1}\bm{v}\mathrm{d}\Omega={\mathcal{N}^{2}\left(\bm{u}\right)}^{\intercal}\left(\int_{\Omega}{\Psi^{2}\left(\bm{x}\right)}^{\intercal}\mathbb{K}^{-1}\left(\bm{x}\right)\Psi^{2}\left(\bm{x}\right)\mathrm{d}\Omega\right)\
\mathcal{N}^{2}\left(\bm{v}\right)={\mathcal{N}^{2}\left(\bm{u}\right)}^{\intercal}\mathbb{M}^{(2)}_{\mathbb{K}^{-1}}\mathcal{N}^{2}\left(\bm{v}\right)\;,$

where
$\mathbb{M}^{(2)}_{\mathbb{K}^{-1}}:=\int_{\Omega}{\Psi^{2}\left(\bm{x}\right)}^{\intercal}\mathbb{K}^{-1}\left(\bm{x}\right)\Psi^{2}\left(\bm{x}\right)\mathrm{d}\Omega$
is the mass matrix of the weighted inner product.

####  3.1.2 Finite dimensional space $S\left(\Omega\right)\subset
L^{2}\left(\Omega\right)$

Let $p$ be any element of the space $S\left(\Omega\right)$ and $N_{v}$ be the
total number of volumes in the discretized domain $\Omega$. If
$\Psi^{3}\left(\bm{x}\right)$ and $\mathcal{N}^{3}\left(p\right)$ form the row
vector of basis functions and the column vector of expansion coefficients,
given as

$\Psi^{3}\left(\bm{x}\right)=\left(\epsilon^{(3)}_{1}\left(\bm{x}\right)\quad\epsilon^{(3)}_{2}\left(\bm{x}\right)\quad\ldots\quad\epsilon^{(3)}_{d_{3}}\left(\bm{x}\right)\right)\;,\qquad\qquad\mathcal{N}^{3}\left(p\right)=\left(\begin{array}[]{c}\mathcal{N}_{1}^{3}\left(p\right)\\\\[4.73611pt]
\mathcal{N}_{2}^{3}\left(p\right)\\\\[4.73611pt] \vdots\\\\[4.73611pt]
\mathcal{N}_{d_{3}}^{3}\left(p\right)\end{array}\right)\;,\qquad\qquad\mbox{respectively,}$

then we can represent $p$ as

$p\left(\bm{x}\right)=\Psi^{3}\left(\bm{x}\right)\mathcal{N}^{3}\left(p\right)\;,$

where $\mathcal{N}^{3}_{i}\left(p\right)=\int_{V_{i}}p\ \mathrm{d}V$ denotes
the integral of $p$ over mesh volume $V_{i}$.

Using this notation, the $L^{2}$-inner product for two elements $p,q\in
S\left(\Omega\right)$ is given as

$\left(p\;,q\right)=\int_{\Omega}p^{\intercal}q\
\mathrm{d}\Omega=\mathcal{N}^{3}\left(p\right)^{\intercal}\left(\int_{\Omega}{\Psi^{3}\left(\bm{x}\right)}^{\intercal}\Psi^{3}\left(\bm{x}\right)\mathrm{d}\Omega\right)\
\mathcal{N}^{3}\left(q\right)=\mathcal{N}^{3}\left(p\right)^{\intercal}\mathbb{M}^{(3)}\
\mathcal{N}^{3}\left(q\right)\;,$ (4)

where,
$\mathbb{M}^{(3)}:=\int_{\Omega}{\Psi^{3}\left(\bm{x}\right)}^{\intercal}\Psi^{3}\left(\bm{x}\right)\mathrm{d}\Omega$
is the mass matrix associated to the basis functions
$\Psi^{3}\left(\bm{x}\right)$.

####  3.1.3 $D_{b}\left(\partial\Omega\right)\subset
H^{-1/2}\left(\partial\Omega\right)$

The trace space $D_{b}\left(\partial\Omega\right)$ is defined as the
restriction of vector fields in $D\left(\Omega\right)$ to the domain boundary
$\partial\Omega$. If $\mathbb{N}_{2}$ is the discrete representation of the
inclusion map, that maps degrees of freedom defined on the boundary to the
global degrees of freedom of the element, and $\mathbb{N}_{2}^{\intercal}$ is
the discrete representation of the trace operator thaat restricts global
degrees of freedom to the boundary, then the basis functions on the boundary
and degrees of freedom on the boundary are given by

$\Psi^{2}_{b}\left(\bm{x}\right):=\Psi^{2}\left(\bm{x}\right)\mathbb{N}_{2}\qquad\mbox{and}\qquad\mathcal{B}^{2}\left(\bm{u}\cdot\bm{n}\right):=\mathbb{N}_{2}^{\intercal}\mathcal{N}^{2}\left(\bm{u}\right)\;.$
(5)

If $\mathrm{tr}\bm{u}\in D_{b}\left(\partial\Omega\right)$ represents the
restriction of flux component of vector field to the domain boundary, then the
solution on the boundary can be represented by

$\mathrm{tr}\
\bm{u}\left(\bm{x}\right)=\Psi^{2}_{b}\left(\bm{x}\right)\mathcal{B}^{2}\left(\bm{u}\right)\;.$

where tr is the trace operator.

For construction of inclusion matrix, see [13, Ex. 3]. The inclusion matrix is
a sparse metric-free matrix, that consists of $+1,-1$ and $0$ entries only,
and is independent of shape and size of elements as long as the topology,
numbering of the degrees of freedom and the orientation of the mesh remains
the same.

####  3.1.4 The divergence operator $\mathbb{E}^{3,2}$

For any element $\bm{u}\in D\left(\Omega\right)$, the divergence operation on
$\bm{u}$ is defined as (see [13, §4.4]):
$\mathrm{div}:D\left(\Omega\right)\longrightarrow S\left(\Omega\right)$, such
that

$\mathrm{div}\
\bm{u}=\Psi^{3}\left(\bm{x}\right)\mathbb{E}^{3,2}\mathcal{N}^{2}\left(\bm{u}\right)\;,$
(6)

where $\mathbb{E}^{3,2}$ is the discrete representation of the divergence
operator that acts on the expansion coefficients of $\bm{u}$. The divergence
operation changes the degrees of freedom and the basis functions to those of
space $S(\Omega)$.

######  Example 1.

Divergence operator for 2D mesh

Figure 1: Left: A square domain divided into $3\times 3$ elements. Right: A
deformed square domain divided into $3\times 3$ elements.

In Figure 1 on the left plot we show a square domain divided into $3\times 3$
elements. The expansion coefficients of $\bm{u}\in D\left(\Omega\right)$, are
defined across the edges in the mesh. The divergence operation on any element
$K$, is then defined as

$\int_{K}\mathrm{div}\ \bm{u}\ \mathrm{d}K=\int_{\partial K}\bm{u}\cdot\bm{n}\
\mathrm{d}\gamma=-\mathcal{N}^{1}\left(\bm{u}\right)_{left}+\mathcal{N}^{1}\left(\bm{u}\right)_{right}-\mathcal{N}^{1}\left(\bm{u}\right)_{bottom}+\mathcal{N}^{1}\left(\bm{u}\right)_{top}\;.$
(7)

If we assemble (7) for all the nine elements, with appropriate numbering, we
get the discrete divergence operator as

$\mathbb{E}^{2,1}=\leavevmode\resizebox{436.8687pt}{}{$\left[\begin{array}[]{cccccccccccccccccccccccc}-1&0&0&1&0&0&0&0&0&0&0&0&-1&0&0&1&0&0&0&0&0&0&0&0\\\
0&0&0&-1&0&0&1&0&0&0&0&0&0&-1&0&0&1&0&0&0&0&0&0&0\\\
0&0&0&0&0&0&-1&0&0&1&0&0&0&0&-1&0&0&1&0&0&0&0&0&0\\\
0&-1&0&0&1&0&0&0&0&0&0&0&0&0&0&-1&0&0&1&0&0&0&0&0\\\
0&0&0&0&-1&0&0&1&0&0&0&0&0&0&0&0&-1&0&0&1&0&0&0&0\\\
0&0&0&0&0&0&0&-1&0&0&1&0&0&0&0&0&0&-1&0&0&1&0&0&0\\\
0&0&-1&0&0&1&0&0&0&0&0&0&0&0&0&0&0&0&-1&0&0&1&0&0\\\
0&0&0&0&0&-1&0&0&1&0&0&0&0&0&0&0&0&0&0&-1&0&0&1&0\\\
0&0&0&0&0&0&0&0&-1&0&0&1&0&0&0&0&0&0&0&0&-1&0&0&1\end{array}\right]\;.$}$

It is a sparse metric free matrix that consists of $+1,-1,0$ entries only. If
the domain is deformed, for example see the right plot of Figure 1, but the
connection between the nodes, edges, surfaces, and volumes, remains the same,
relation (7) remains the same and consequently the matrix $\mathbb{E}^{2,1}$
remains unchanged.

###  3.2 Algebraic dual representations

In this section we will introduce the finite dimensional sub spaces for
$L^{2}\left(\Omega\right)$ and $H^{1/2}\left(\partial\Omega\right)$ using the
algebraic dual representations as defined in [13].

####  3.2.1 Finite dimensional space $\widetilde{S}\left(\Omega\right)\subset
L^{2}\left(\Omega\right)$

For the finite dimensional space $S\left(\Omega\right)$ let
$\widetilde{S}\left(\Omega\right)$ be the corresponding algebraic dual
representation . For any element
${\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}q}\in\widetilde{S}\left(\Omega\right)$,
if $\widetilde{\Psi}^{0}\left(\bm{x}\right)$ and
$\widetilde{\mathcal{N}}^{0}\left(\widetilde{q}\right)$ are the associated
basis functions and the expansion coefficients, then we can represent $q$ as

${\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}q}\left(\bm{x}\right)=\widetilde{\Psi}^{0}\left(\bm{x}\right)\widetilde{\mathcal{N}}^{0}\left({\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}q}\right)\;,$

where

$\widetilde{\Psi}^{0}\left(\bm{x}\right)=\Psi^{3}\left(\bm{x}\right)\left(\mathbb{M}^{(3)}\right)^{-1}\;,\qquad\text{and}\;,\qquad\widetilde{\mathcal{N}}^{0}\left({\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}q}\right)=\mathbb{M}^{(3)}\mathcal{N}^{3}\left(q\right)\;.$

Using Corollary 1, the $L^{2}$-inner product between the elements $p,q\in
S\left(\Omega\right)\times\widetilde{S}\left(\Omega\right)$ is then given by

$\left(q\;,p\right)=\int_{\Omega}q^{\intercal}p\
\mathrm{d}\Omega={\widetilde{\mathcal{N}}^{0}\left({\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}q}\right)}^{\intercal}\mathcal{N}^{3}\left(p\right)\;.$
(8)

We see that the $L^{2}$-inner product in (4) requires the evaluation of mass
matrix, whereas the inner product in (8) requires only the vector product
between the expansion coefficients which makes the discrete system more sparse
and easier to set up.

####  3.2.2 Finite dimensional trace space
$\widetilde{D}_{b}\left(\partial\Omega\right)\subset
H^{1/2}\left(\partial\Omega\right)$

It is known that the Sobolev spaces $H^{1/2}\left(\partial\Omega\right)$ and
$H^{-1/2}\left(\partial\Omega\right)$ are dual to each other. To replicate
this duality also in the discrete setting, we choose
$\widetilde{D}_{b}\left(\partial\Omega\right)$ as finite dimensional sub space
of $H^{1/2}\left(\partial\Omega\right)$ which is the algebraic dual
representation of $D_{b}\left(\partial\Omega\right)\subset
H^{-1/2}\left(\partial\Omega\right)$.

For an element $\lambda\in\widetilde{D}_{b}\left(\partial\Omega\right)$ if
$\widetilde{\Psi}_{b}^{0}\left(\bm{x}\right)$ and
$\widetilde{\mathcal{B}}^{0}\left(\lambda\right)$ are the basis functions and
the expansion coefficients then, we can represent $\lambda$ as

$\lambda\left(\bm{x}\right)=\widetilde{\Psi}^{0}_{b}\left(\bm{x}\right)\widetilde{\mathcal{B}}^{0}\left(\lambda\right)\;.$
(9)

where

$\widetilde{\Psi}^{0}_{b}\left(\bm{x}\right)=\Psi^{2}_{b}\left(\bm{x}\right)\left(\mathbb{M}_{b}^{(2)}\right)^{-1}\qquad\widetilde{\mathcal{B}}^{0}\left(\lambda\right)=\mathbb{M}^{(2)}_{b}\mathcal{B}^{2}\left(\lambda\right)\;,\qquad\mbox{with,}\qquad\mathbb{M}_{b}^{(2)}=\int\Psi^{2}_{b}\left(\bm{x}\right)^{\intercal}\Psi^{2}_{b}\left(\bm{x}\right)\mathrm{d}\Gamma\;.$
(10)

The inner product between the elements $\lambda\;,\mathrm{tr}\
\bm{u}\in\widetilde{D}_{b}\left(\partial\Omega\right)\times
D_{b}\left(\partial\Omega\right)$, is given by

$\left(\lambda\;,\mathrm{tr}\
\bm{u}\right)=\widetilde{\mathcal{B}}^{0}\left(\lambda\right)^{\intercal}\left(\int_{\Omega}\widetilde{\Psi}^{0}_{b}\left(\bm{x}\right)^{\intercal}\Psi^{2}_{b}\left(\bm{x}\right)\mathrm{d}\Omega\right)\mathcal{B}^{2}\left(\mathrm{tr}\
\bm{u}\right)=\widetilde{\mathcal{B}}^{0}\left(\lambda\right)^{\intercal}\mathcal{B}^{2}\left(\mathrm{tr}\
\bm{u}\right)\;.$ (11)

###  3.3 The gradient of dual representations

For a scalar field $p\in\widetilde{S}\left(\Omega\right),$ and
$\hat{p}\in\widetilde{D}_{b}\left(\partial\Omega\right)$ its values on the
domain boundary, the gradient operation for dual representations is defined as
[13, Def. 19],
$\widetilde{\mathrm{grad}}:\widetilde{S}\left(\Omega\right)\times\widetilde{D}_{b}\left(\partial\Omega\right)\rightarrow\widetilde{D}\left(\Omega\right)$,
such that

$\int_{\Omega}\widetilde{\mathrm{grad}}\left(p\;,\hat{p}\right)\bm{u}\
\mathrm{d}\Omega=-\int_{\Omega}p\left(\mathrm{div}\ \bm{u}\right)\
\mathrm{d}\Omega+\int_{\partial\Omega}\hat{p}\left(\bm{u}\cdot\bm{n}\right)\
\mathrm{d}\Gamma\qquad\qquad\forall\bm{u}\in D\left(\Omega\right)\;.$ (12)

The expansion coefficients of $\widetilde{\mathrm{grad}}\
\left(p,\hat{p}\right)$ are given as

$\widetilde{\mathcal{N}}^{1}\left(\widetilde{\mathrm{grad}}\left({p}\;,\hat{p}\right)\right)=-{\mathbb{E}^{3,2}}^{\intercal}\widetilde{\mathcal{N}}^{0}\left(p\right)+\mathbb{N}_{2}\widetilde{\mathcal{B}}^{0}\left(\hat{p}\right)\;,$
(13)

expanded in the basis
$\widetilde{\Psi}^{1}\left(\bm{x}\right)=\mathbb{M}^{(2)}\Psi^{2}\left(\bm{x}\right)$.

###  3.4 Broken finite dimensional spaces

In this section we define the finite dimensional subset of Broken Sobolev
spaces.

####  3.4.1 Finite dimensional space $D\left(\Omega_{T}\right)\subset
H\left(\mathrm{div};\Omega_{T}\right)$

For the set of domains $\Omega_{T}$ we define the finite element space
$D\left(\Omega_{T}\right)\subset H\left(\mathrm{div};\Omega_{T}\right)$ as

$D\left(\Omega_{T}\right):=\prod_{\mathcal{M}\in\Omega_{T}}D\left(\mathcal{M}\right)\;.$
(14)

For any two elements $\bm{u}\;,\bm{v}\in D\left(\Omega_{T}\right)$, and
symmetric positive definite permeability tensor $\mathbb{K}$, the weighted
$L^{2}$-inner product is given by

$\left(\bm{u}\;,\bm{v}\right)_{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\mathbb{K}^{-1}}}=\int_{\Omega}\bm{u}^{\intercal}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\mathbb{K}^{-1}}\bm{v}\
\mathrm{d}\Omega=\sum_{i=1}^{T}\int_{\mathcal{M}_{i}}\bm{u}^{\intercal}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\mathbb{K}^{-1}}\bm{v}\
\mathrm{d}\mathcal{M}={{\mathcal{N}}^{2}\left(\bm{u}\right)}^{\intercal}{\mathbb{M}}^{(2)}_{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\mathbb{K}^{-1}}}{\mathcal{N}}^{2}\left(\bm{v}\right)\;,$
(15)

where, $\mathcal{N}^{2}\left(\bm{u}\right)$ and
$\mathcal{N}^{2}\left(\bm{v}\right)$ are the column vector of assembled
expansion coefficients of all the sub domains and the mass matrix
$\mathbb{M}^{(2)}_{\mathbb{K}^{-1}}$ is given as

${\mathbb{M}}^{(2)}_{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\mathbb{K}^{-1}}}=\left[\begin{array}[]{ccccc}\mathbb{M}^{(2)}_{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\mathbb{K}^{-1}},1}\\\\[4.73611pt]
&\mathbb{M}^{(2)}_{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\mathbb{K}^{-1}},2}\\\\[4.73611pt]
&&\ddots\\\\[4.73611pt]
&&&\mathbb{M}^{(2)}_{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\mathbb{K}^{-1}},T}\end{array}\right]\;,\qquad\text{with,}\qquad\mathbb{M}^{(2)}_{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\mathbb{K}^{-1}},i}=\int_{\mathcal{M}_{i}}{\Psi^{2}_{i}\left(\bm{x}\right)}^{\intercal}{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\mathbb{K}^{-1}}\left(\bm{x}\right)}\Psi^{2}_{i}\left(\bm{x}\right)\mathrm{d}\mathcal{M}\;,$
(16)

where $\Psi_{i}^{2}\left(\bm{x}\right)$ are the basis functions associated to
the finite dimensional space $D\left(\mathcal{M}_{i}\right)$.

####  3.4.2 Finite dimensional space $S\left(\Omega_{T}\right)\subset
L^{2}\left(\Omega_{T}\right)$

For the set of domains $\Omega_{T}$ we define the finite element space
$S\left(\Omega_{T}\right)\subset L^{2}\left(\Omega_{T}\right)$ as

$S\left(\Omega_{T}\right):=\prod_{\mathcal{M}\in\Omega_{T}}S\left(\mathcal{M}\right)\;.$
(17)

For elements $p,q\in S\left(\Omega_{T}\right)$, the $L^{2}$\- inner product is
given by

$\left(p,q\right)=\int_{\Omega}p^{\intercal}q\
\mathrm{d}\Omega=\sum_{i=1}^{T}\int_{\mathcal{M}_{i}}p^{\intercal}q\
\mathrm{d}\mathcal{M}={\mathcal{N}^{3}\left(p\right)}^{\intercal}\mathbb{M}^{(3)}\mathcal{N}^{3}\left(q\right)\;,$
(18)

where $\mathcal{N}^{3}\left(p\right)$, $\mathcal{N}^{3}\left(q\right)$ are the
column vectors of assembled expansion coefficients of all the sub domains, and
the mass matrix $\mathbb{M}^{(3)}$ is the block diagonal mass matrix given by

$\mathbb{M}^{(3)}=\left[\begin{array}[]{ccccc}\mathbb{M}^{(3)}_{1}\\\\[4.73611pt]
&\mathbb{M}^{(3)}_{2}\\\\[4.73611pt] &&\ddots\\\\[4.73611pt]
&&&\mathbb{M}^{(3)}_{T}\end{array}\right]\;,\qquad\text{with,}\qquad\mathbb{M}^{(3)}_{i}=\int_{\mathcal{M}_{i}}{\Psi^{3}_{i}\left(\bm{x}\right)}^{\intercal}\Psi^{3}_{i}\left(\bm{x}\right)\mathrm{d}\mathcal{M}\;,$
(19)

where $\Psi_{i}^{3}\left(\bm{x}\right)$ are the basis associated with the
finite dimensional space $S\left(\mathcal{M}_{i}\right)$.

####  3.4.3 Finite dimensional space
$D_{b}\left(\partial\Omega_{T}\right)\subset
H^{-1/2}\left(\partial\Omega_{T}\right)$

For the set of boundaries $\partial\Omega_{T}$ we define the trace space
$D_{b}\left(\partial\Omega_{T}\right)\subset
H^{-1/2}\left(\partial\Omega_{T}\right)$ as

$D_{b}\left(\partial\Omega_{T}\right):=\prod_{\mathcal{M}\in\Omega_{T}}D_{b}\left(\partial\mathcal{M}\right)\;.$
(20)

####  3.4.4 Finite dimensional space
$\widetilde{S}\left(\Omega_{T}\right)\subset L^{2}\left(\Omega_{T}\right)$

The algebraic dual representation of $S\left(\Omega_{T}\right)$ is denoted by
$\widetilde{S}\left(\Omega_{T}\right)$ and defined as

$\widetilde{S}\left(\Omega_{T}\right):=\prod_{\mathcal{M}\in\Omega_{T}}\widetilde{S}\left(\mathcal{M}\right)\;.$
(21)

Let $p,q\in
S\left(\Omega_{T}\right)\times\widetilde{S}\left(\Omega_{T}\right)$, then the
inner product between the elements is given by

$\left(p,q\right)=\sum_{i=1}^{T}\int_{\mathcal{M}_{i}}q^{\intercal}p\
\mathrm{d}\mathcal{M}={\widetilde{\mathcal{N}}^{0}\left(q\right)}^{\intercal}\mathcal{N}^{3}\left(p\right)\;,$
(22)

where,
$\widetilde{\mathcal{N}}^{0}\left(q\right)\;,\mathcal{N}^{3}\left(p\right)$
are the column vectors of assembled expansion coefficients.

####  3.4.5 Finite dimensional space
$\widetilde{D}_{b}\left(\partial\Omega_{e}\right)\subset
H^{1/2}\left(\partial\Omega_{e}\right)$

Let $\gamma\in\partial\Omega_{e}$ be any boundary face of the sub domains. Let
$D_{b}\left(\gamma\right)$ be the space of restriction of basis of
$D\left(\Omega_{T}\right)$ to $\gamma$ and
$\widetilde{D}_{b}\left(\gamma\right)$ be its algebraic dual representation .

Then the finite dimensional sub space
$\widetilde{D}_{b}\left(\partial\Omega_{e}\right)\subset
H^{1/2}\left(\partial\Omega_{e}\right)$ is defined as

$\widetilde{D}_{b}\left(\partial\Omega_{e}\right)=\prod_{\gamma\in\partial\Omega_{e}}\widetilde{D}_{b}\left(\gamma\right)\;.$
(23)

For elements $\lambda,\mathrm{tr}\
\bm{u}\in\widetilde{D}_{b}\left(\partial\Omega_{e}\right)\times
D_{b}\left(\partial\Omega_{T}\right)$, we have the inner product given by

$\left(\lambda,\mathrm{tr}\
\bm{u}\right)=\sum_{i=1}^{T}\int_{\partial\mathcal{M}_{i}}\lambda^{\intercal}\mathrm{tr}\
\bm{u}\
\mathrm{d}\Gamma=\sum_{i=1}^{T}\widetilde{\mathcal{B}}_{i}^{0}\left(\lambda\right)^{\intercal}\mathcal{B}_{i}^{2}\left(\mu\right)\;,$
(24)

where $\widetilde{\mathcal{B}}^{0}\left(\lambda\right)$,
$\mathcal{B}^{2}\left(\mu\right)$ are the assembled coefficients over all sub
domains.

If $\mu=\bm{u}\cdot\bm{n}|_{\partial\mathcal{M}_{i}}$ for $i=1,\ldots,T$, i.e.
flux component at sub domain boundaries, then using (5) we have that

$\mathcal{B}^{2}\left(\mu\right)=\mathcal{B}^{2}\left(\bm{u}\cdot\bm{n}\right)=\mathbb{N}_{2}^{\intercal}\mathcal{N}^{2}\left(\bm{u}\right)\;,$
(25)

and we can write (24) as

$\left(\lambda,\bm{u}\cdot\bm{n}\right)=\widetilde{B}^{0}\left(\lambda\right)^{\intercal}\mathbb{N}_{2}^{\intercal}\mathcal{N}^{2}\left(\bm{u}\right)\;,$
(26)

where, $\mathcal{N}^{2}\left(\bm{u}\right)$ is the column vector of assembled
expansion coefficients and $\mathbb{N}_{2}^{\intercal}$ is the assembled trace
matrix of all the sub domains.

####  3.4.6 The divergence operator for broken Sobolev spaces

Let $q,\bm{u}\in\widetilde{S}\left(\Omega_{T}\right)\times
D\left(\Omega_{T}\right)$ be the elements of the broken finite dimensional
spaces. Then the divergence operation on the vector field $\bm{u}$ is given by

$\left(q\;,\mathrm{div}\ \bm{u}\right)=\sum_{i=1}^{T}\int_{\mathcal{M}_{i}}q\
\mathrm{div}\ \bm{u}\
\mathrm{d}\mathcal{M}={\widetilde{\mathcal{N}}^{0}\left(q\right)}^{\intercal}\mathbb{E}^{3,2}\mathcal{N}^{2}\left(\bm{u}\right)\qquad\forall
q\in\widetilde{S}\left(\Omega_{T}\right)\;,$ (27)

where,
$\widetilde{\mathcal{N}}^{0}\left(q\right),\mathcal{N}^{2}\left(\bm{u}\right)$
are the column vectors of the assembled expansion coefficients, and
$\mathbb{E}^{3,2}$ is the assembled divergence operator. If
$\mathbb{E}^{3,2}_{i}$ is the discrete representation of divergence operator
on the domain $\mathcal{M}_{i}$, then we have

$\mathbb{E}^{3,2}=\left[\begin{array}[]{ccccc}\mathbb{E}^{3,2}_{1}\\\\[4.73611pt]
&\mathbb{E}^{3,2}_{2}\\\\[4.73611pt] &&\ddots\\\\[4.73611pt]
&&&\mathbb{E}^{3,2}_{T}\end{array}\right]\;.$ (28)

If the topology of sub domain discretizations is also the same, see Example 1,
then we have that
$\mathbb{E}^{3,2}_{1}=\mathbb{E}^{3,2}_{2}=\ldots=\mathbb{E}^{3,2}_{T}$. In
this paper we have only used the case with the same topology for all the sub
domains, therefore

$\mathbb{E}^{3,2}=\left[\begin{array}[]{ccccc}\mathbb{E}^{3,2}_{i}\\\\[4.73611pt]
&\mathbb{E}^{3,2}_{i}\\\\[4.73611pt] &&\ddots\\\\[4.73611pt]
&&&\mathbb{E}^{3,2}_{i}\end{array}\right]\;,$ (29)

where $\mathbb{E}^{3,2}_{i}$ is the topological divergence operator for any of
the sub domains.

####  3.4.7 The gradient of dual representations

Let $p$ be a scalar field represented by algebraic dual representation
$\widetilde{S}\left(\Omega_{T}\right)$, and $\hat{p}$ be the boundary value of
the scalar field on the sub domain boundary faces
$\gamma\in\partial\Omega_{e}$. The gradient operation for dual representations
$p\;,\hat{p}\in\widetilde{S}\left(\Omega_{T}\right)\times\widetilde{D}_{b}\left(\partial\Omega_{e}\right)$
is then defined as,
$\widetilde{\mathrm{grad}}:\widetilde{S}\left(\Omega_{T}\right)\times\widetilde{D}_{b}\left(\partial\Omega_{e}\right)\rightarrow\widetilde{D}\left(\Omega_{T}\right)$,
such that

$\int_{\Omega}\widetilde{\mathrm{grad}}\left(p\;,\hat{p}\right)\bm{u}\
\mathrm{d}\Omega=\sum_{i=1}^{T}\int_{\mathcal{M}_{i}}\widetilde{\mathrm{grad}}\left(p\;,\hat{p}\right)\bm{u}\
\mathrm{d}\mathcal{M}=\sum_{i=1}^{T}\left(-\int_{\mathcal{M}_{i}}p\left(\mathrm{div}\
\bm{u}\right)+\int_{\partial\mathcal{M}_{i}}\hat{p}\left(\bm{u}\cdot\bm{n}\right)\
\mathrm{d}\Gamma\right)\quad\forall\bm{u}\in D\left(\Omega_{T}\right)\;.$ (30)

For the pair of elements
$p\;,\hat{p}\in\widetilde{S}\left(\Omega_{T}\right)\times\widetilde{D}\left(\partial\Omega_{e}\right)$
the $H^{1}$-norm is then defined as

$\|p\|^{2}_{H^{1}\left(\Omega\right)}=\|p\|^{2}_{L^{2}\left(\Omega\right)}+\|\widetilde{\mathrm{grad}}\left(p,\hat{p}\right)\|^{2}_{L^{2}\left(\Omega\right)}=\sum_{i=1}^{T}\left\\{\|p\|^{2}_{L^{2}\left(\mathcal{M}_{i}\right)}+\|\widetilde{\mathrm{grad}}\left(p,\hat{p}\right)\|^{2}_{L^{2}\left(\mathcal{M}_{i}\right)}\right\\}\;.$
(31)

##  4 Model anisotropic diffusion problem

In this section we will use the broken finite dimensional spaces defined in
Section 3 to derive the algebraic formulation for DD formulation of Darcy
problem.

The equations for Darcy problem in the domain $\Omega$ are given by

$\left\\{\begin{array}[]{ll}\bm{u}+\mathbb{K}\ \mathrm{grad}\
p&=0\\\\[4.73611pt] \mathrm{div}\
\bm{u}&=f\end{array}\right.\quad\mbox{with}\quad\left\\{\begin{array}[]{ll}\bm{u}\cdot\bm{n}=\hat{u}&\mbox{on}\quad\Gamma_{N}\\\\[4.73611pt]
p=\hat{p}&\mbox{on}\quad\Gamma_{D}\\\\[4.73611pt]
\Gamma_{N}\cap\Gamma_{D}=\emptyset&\mbox{and}\quad\Gamma_{N}\cup\Gamma_{D}=\partial\Omega\end{array}\right.\;,$
(32)

where $\bm{u}$ is the velocity, $\mathbb{K}$ is the symmetric positive
definite permeability tensor, $p$ is the pressure, $f$ is the given right hand
side term, $\bm{n}$ is the outward unit normal vector, $\hat{u}$ is the given
velocity boundary condition imposed on Neumann boundary $\Gamma_{N}$ and
$\hat{p}$ is the pressure boundary condition imposed on the Dirichlet boundary
$\Gamma_{D}$.

The Lagrange functional for continuous formulation of Darcy equations is given
by

$\mathcal{L}\left(\bm{u},p;f,\hat{p},\hat{u}\right)=\int_{\Omega}\frac{1}{2}\bm{u}^{\intercal}\mathbb{K}^{-1}\bm{u}\
\mathrm{d}\Omega-\int_{\Omega}p\left(\mathrm{div}\
\bm{u}-f\right)\mathrm{d}\Omega+\int_{\Gamma_{D}}\hat{p}\left(\bm{u}\cdot\bm{n}\right)\mathrm{d}\Gamma+\int_{\Gamma_{N}}\left(\bm{u}\cdot\bm{n}-\hat{u}\right)\mathrm{d}\Gamma\;.$
(33)

The algebraic system for (33), using algebraic dual representations of [13],
is given by

$\left[\begin{array}[]{cc}\mathbb{M}_{\mathbb{K}^{-1}}^{(2)}&-{\mathbb{E}^{3,2}}^{\intercal}\\\\[4.73611pt]
-\mathbb{E}^{3,2}&0\end{array}\right]\left[\begin{array}[]{c}\mathcal{N}^{2}\left(\bm{u}\right)\\\\[4.73611pt]
\widetilde{\mathcal{N}}^{0}\left(p\right)\end{array}\right]=\left[\begin{array}[]{c}-{\mathbb{N}_{2}}\widetilde{\mathcal{B}}^{0}\left(\hat{p}\right)\\\\[4.73611pt]
-\mathcal{N}^{3}\left(f\right)\end{array}\right]\;.$ (34)

The above system can be solved for unknowns of $p$, and the unknowns of
$\bm{u}$ as

$\displaystyle\widetilde{\mathcal{N}}^{0}\left(p\right)$ $\displaystyle=$
$\displaystyle\left(\mathbb{E}^{3,2}{\mathbb{M}^{(2)}_{\mathbb{K}^{-1}}}^{-1}{\mathbb{E}^{3,2}}^{\intercal}\right)^{-1}\left(\mathcal{N}^{3}\left(f\right)+\mathbb{E}^{3,2}{\mathbb{M}^{(2)}_{\mathbb{K}^{-1}}}^{-1}\mathbb{N}_{2}\widetilde{\mathcal{B}}^{0}\left(\hat{p}\right)\right)$
(35) $\displaystyle\mathcal{N}^{2}\left(\bm{u}\right)$ $\displaystyle=$
$\displaystyle{\mathbb{M}_{\mathbb{K}^{-1}}^{(2)}}^{-1}\left({\mathbb{E}^{3,2}}^{\intercal}\widetilde{\mathcal{N}}^{0}\left(p\right)-\mathbb{N}_{2}\widetilde{\mathcal{B}}^{0}\left(\hat{p}\right)\right)$
(36)

For the DD formulation of (32) we break the domain $\Omega$ into $T$ sub
domains, see (1), and we use the Lagrange multipliers $\lambda$ to enforce the
required continuity across the sub domains. The weak formulation for (32) is
then obtained using the Lagrange functional

$\displaystyle\mathcal{L}\left(\bm{u},p,\lambda;f,\hat{p},\hat{u}\right)$
$\displaystyle=$
$\displaystyle\sum_{i=1}^{T}\left\\{\int_{\mathcal{M}_{i}}\frac{1}{2}\bm{u}^{\intercal}\mathbb{K}^{-1}\bm{u}\
\mathrm{d}\Omega-\int_{\mathcal{M}_{i}}p\left(\mathrm{div}\
\bm{u}-f\right)\mathrm{d}\Omega+\int_{\partial\mathcal{M}_{i}\setminus(\Gamma_{D}\cup\Gamma_{N})}\lambda\left(\bm{u}\cdot\bm{n}\right)\mathrm{d}\Gamma\right.$
(37)
$\displaystyle\left.+\int_{\partial\mathcal{M}_{i}\cap\Gamma_{N}}\lambda\left(\bm{u}\cdot\bm{n}-\hat{u}\right)\mathrm{d}\Gamma+\int_{\partial\mathcal{M}_{i}\cap\Gamma_{D}}\hat{p}\left(\bm{u}\cdot\bm{n}\right)\mathrm{d}\Gamma\right\\}\;,$

where, on the left hand side, the first term is the kinetic energy term, the
second term imposes the constraint on divergence of velocity field $\bm{u}$,
the third term imposes continuity of flux across the sub domain faces, the
fourth term imposes the Neumann boundary condition and the fifth term imposes
the Dirichlet boundary conditions.

######  Lemma 2.

The Lagrange multipliers are pressure boundary values on sub domains.

###### Proof.

In (37) if we take variations with respect to $\bm{u}$, we get

$\sum_{i=1}^{T}\left\\{\int_{\mathcal{M}_{i}}\bm{v}^{\intercal}\mathbb{K}^{-1}\bm{u}\
\mathrm{d}\mathcal{M}-\int_{\mathcal{M}_{i}}\left(\mathrm{div}\
\bm{v}\right)p\
\mathrm{d}\mathcal{M}+\int_{\partial\mathcal{M}_{i}\setminus{\Gamma_{D}}}\left(\bm{v}\cdot\bm{n}\right)\lambda\
\mathrm{d}\Gamma+\int_{\partial\mathcal{M}_{i}\cap\Gamma_{D}}\hat{p}\left(\bm{v}\cdot\bm{n}\right)\mathrm{d}\Gamma\right\\}=0\qquad\forall\bm{v}\in
H\left(\mathrm{div};\Omega_{t}\right)\;.$

For any sub domain $\mathcal{M}_{i}$, we have

$\int_{\mathcal{M}_{i}}\bm{v}^{\intercal}\mathbb{K}^{-1}\bm{u}\
\mathrm{d}\mathcal{M}-\int_{\mathcal{M}_{i}}\left(\mathrm{div}\
\bm{v}\right)p\
\mathrm{d}\mathcal{M}+\int_{\partial\mathcal{M}_{i}\setminus{\Gamma_{D}}}\left(\bm{v}\cdot\bm{n}\right)\lambda\
\mathrm{d}\Gamma+\int_{\partial\mathcal{M}_{i}\cap\Gamma_{D}}\left(\bm{v}\cdot\bm{n}\right)\hat{p}\mathrm{d}\Gamma=0\qquad\forall\bm{v}\in
H\left(\mathrm{div};\mathcal{M}_{i}\right)\;.$

If the solution is sufficiently smooth, then using integration by parts on the
second term, we get

$\int_{\mathcal{M}_{i}}\bm{v}^{\intercal}\mathbb{K}^{-1}\bm{u}\
\mathrm{d}\mathcal{M}-\left(-\int_{\mathcal{M}_{i}}\bm{v}^{\intercal}\left(\mathrm{grad}\
p\right)\mathrm{d}\mathcal{M}+\int_{\partial\mathcal{M}_{i}}\left(\bm{v}\cdot\bm{n}\right)p\
\mathrm{d}\Gamma\right)+\int_{\partial\mathcal{M}_{i}\setminus{\Gamma_{D}}}\left(\bm{v}\cdot\bm{n}\right)\lambda\
\mathrm{d}\Gamma+\int_{\partial\mathcal{M}_{i}\cap\Gamma_{D}}\left(\bm{v}\cdot\bm{n}\right)\hat{p}\
\mathrm{d}\Gamma=0\quad\forall\bm{v}\in
H\left(\mathrm{div};\mathcal{M}_{i}\right)\;.$

Now, combining the first and the second term, and the third, fourth and fifth
term we get

$\int_{\mathcal{M}_{i}}\bm{v}^{\intercal}\left(\mathbb{K}^{-1}\bm{u}+\mathrm{grad}\
p\right)\
\mathrm{d}\mathcal{M}+\int_{\partial\mathcal{M}_{i}\setminus{\Gamma_{D}}}\left(\bm{v}\cdot\bm{n}\right)\left(\lambda-p\right)\mathrm{d}\Gamma+\int_{\partial\mathcal{M}_{i}\cap\Gamma_{D}}\left(\bm{v}\cdot\bm{n}\right)\left(p-\hat{p}\right)\mathrm{d}\Gamma=0\qquad\forall\bm{v}\in
H\left(\mathrm{div};\mathcal{M}_{i}\right)\;.$

This implies that $\mathbb{K}^{-1}\bm{u}+\mathrm{grad}\ p=0$ in the
$L^{2}$-sense, i.e. almost everywhere, and $\lambda-p=0$ in the
$H^{1/2}$-sense, i.e. between the sub domains and $p=\hat{p}$ along
$\Gamma_{D}$. As this should hold for all $\bm{v}\in
H\left(\mathrm{div};\mathcal{M}_{i}\right)$, we have that $\lambda=p$, i.e.
the Lagrange multipliers are the pressure boundary values on the sub domain
boundaries, $\partial\mathcal{M}_{i}$. ∎

The optimality conditions at continuous level for the Lagrange functional (37)
are given by: For given $f\in L^{2}\left(\Omega\right)$, $\hat{p}\in
H^{1/2}\left(\Gamma_{D}\right)$, $\hat{u}\in H^{-1/2}\left(\Gamma_{N}\right)$,
find $\bm{u}\in H\left(\mathrm{div};\mathcal{M}_{i}\right)$, $p\in
L^{2}\left(\mathcal{M}_{i}\right)$, $\lambda\in
H^{1/2}\left(\partial\Omega_{e}\right)$, such that

$\sum_{i=1}^{T}\left\\{\begin{array}[]{clllll}\left(\bm{v},\mathbb{K}^{-1}\bm{u}\right)&-\left(\mathrm{div}\
\bm{v},p\right)&+\left(\bm{v}\cdot\bm{n},\lambda\right)&=-\left(\bm{v}\cdot\bm{n},\hat{p}\right)&\forall\bm{v}\in
H\left(\mathrm{div};\mathcal{M}_{i}\right)\\\\[4.73611pt]
-\left(q,\mathrm{div}\ \bm{u}\right)&&&=-\left(q,f\right)&\forall q\in
L^{2}\left(\mathcal{M}_{i}\right)\\\\[4.73611pt]
\left(\mu,\bm{u}\cdot\bm{n}\right)&&&=\left(\mu,\hat{u}\right)&\forall\mu\in
H^{1/2}\left(\partial\Omega_{e}\right)\end{array}\right.\;.$ (38)

The finite dimensional problem is then given by: For given $f\in
L^{2}\left(\Omega\right)$, $\hat{p}\in H^{1/2}\left(\Gamma_{D}\right)$,
$\hat{u}\in H^{-1/2}\left(\Gamma_{N}\right)$, find $\bm{u}\in
D\left(\mathcal{M}_{i}\right)$,
${p}\in\widetilde{S}\left(\mathcal{M}_{i}\right)$,
${\lambda}\in\widetilde{D}_{b}\left(\partial\Omega_{e}\right)$, such that

$\sum_{i=1}^{T}\left\\{\begin{array}[]{clllll}\left(\bm{v},\mathbb{K}^{-1}\bm{u}\right)&-\left(\mathrm{div}\
\bm{v},{p}\right)&+\left(\bm{v}\cdot\bm{n},{\lambda}\right)&=-\left(\bm{v}\cdot\bm{n},{\hat{p}}\right)&\forall\bm{v}\in
D\left(\mathcal{M}_{i}\right)\\\\[4.73611pt] -\left({q},\mathrm{div}\
\bm{u}\right)&&&=-\left({q},f\right)&\forall
q\in\widetilde{S}\left(\mathcal{M}_{i}\right)\\\\[4.73611pt]
\left({\mu},\bm{u}\cdot\bm{n}\right)&&&=\left({\mu},\hat{u}\right)&\forall{\mu}\in\widetilde{D}_{b}\left(\partial\Omega_{e}\right)\end{array}\right.\;.$
(39)

In (39) we see that all terms, except one - the weighted inner product term,
are the inner product between the primal and the dual representations.
Consequently, these terms do not require evaluation of (dense) mass matrices.
The only matrices associated with these terms are the sparse, metric-free
divergence operator, or the inclusion matrix. See also (40).

###  4.1 The algebraic formulation

The inner product for the left hand side terms of (39) are evaluated as

$\displaystyle\left(\bm{v},\mathbb{K}^{-1}\bm{u}\right)$
$\displaystyle\stackrel{{\scriptstyle\eqref{eq:inner_D}}}{{=}}$
$\displaystyle{\mathcal{N}^{2}\left(\bm{v}\right)}^{\intercal}\mathbb{M}^{(2)}_{\mathbb{K}^{-1}}\mathcal{N}^{2}\left(\bm{u}\right)\;,$
$\displaystyle\left({q}\;,\mathrm{div}\ \bm{u}\right)$
$\displaystyle\stackrel{{\scriptstyle\eqref{eq:inner_div}}}{{=}}$
$\displaystyle{\widetilde{\mathcal{N}}^{0}\left({q}\right)}^{\intercal}\mathbb{E}^{3,2}\mathcal{N}^{2}\left(\bm{u}\right)\;,$
$\displaystyle\left({\mu}\;,\bm{u}\cdot\bm{n}\right)$
$\displaystyle\stackrel{{\scriptstyle\eqref{eq:inner_trace}}}{{=}}$
$\displaystyle{\widetilde{\mathcal{B}}^{0}\left({\mu}\right)}^{\intercal}\mathbb{N}_{2}^{\intercal}\mathcal{N}^{2}\left(\bm{u}\right)\;,$

and the inner product for the right hand side terms of (39) are evaluated as

$\displaystyle\left(\bm{v}\cdot\bm{n}\;,{\hat{p}}\right)$
$\displaystyle\stackrel{{\scriptstyle\eqref{eq:inner_trace}}}{{=}}$
$\displaystyle{\mathcal{N}^{2}\left(\bm{v}\right)}^{\intercal}\mathbb{N}_{2}\widetilde{\mathcal{B}}^{0}\left({\hat{p}}\right){\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}={\mathcal{N}^{2}\left(\bm{v}\right)}^{\intercal}\mathbb{N}_{2}\int_{\Gamma_{D}}\widetilde{\Psi}^{0}_{b}\left(\bm{x}\right)\hat{p}\left(\bm{x}\right)\mathrm{d}\Gamma}\;,$
$\displaystyle\left({q}\;,f\right)$
$\displaystyle\stackrel{{\scriptstyle\eqref{eq:inner_S_dual}}}{{=}}$
$\displaystyle{\widetilde{\mathcal{N}}^{0}\left({q}\right)}^{\intercal}\mathcal{N}^{3}\left(f\right){\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}={\widetilde{\mathcal{N}}^{0}\left({q}\right)}^{\intercal}\int_{\Omega}\widetilde{\Psi}^{0}\left(\bm{x}\right)f\left(x\right)\mathrm{d}\Omega}\;,$
$\displaystyle\left({\mu}\;,\hat{u}\right)$
$\displaystyle\stackrel{{\scriptstyle\eqref{eq:inner_trace2}}}{{=}}$
$\displaystyle{\widetilde{\mathcal{B}}^{0}\left({\mu}\right)}^{\intercal}\mathcal{B}^{2}\left(\hat{u}\right){\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}={\widetilde{\mathcal{B}}^{0}\left({\mu}\right)}^{\intercal}\int_{\Gamma_{N}}\Psi_{b}^{2}\left(\bm{x}\right)\hat{u}\left(\bm{x}\right)\mathrm{d}\Gamma}\;.$

Using above relations we can write the algebraic formulation for (39) as

$\left[\begin{array}[]{lcc}\mathbb{M}^{(2)}_{\mathbb{K}^{-1}}&-{\mathbb{E}^{3,2}}^{\intercal}&\mathbb{N}_{2}\\\\[4.73611pt]
-\mathbb{E}^{3,2}&0&0\\\\[4.73611pt]
\mathbb{N}_{2}^{\intercal}&0&0\end{array}\right]\left[\begin{array}[]{c}\mathcal{N}^{2}\left(\bm{u}\right)\\\\[4.73611pt]
\widetilde{\mathcal{N}}^{0}\left({p}\right)\\\\[4.73611pt]
\widetilde{\mathcal{B}}^{0}\left({\lambda}\right)\end{array}\right]=\left[\begin{array}[]{c}-\mathbb{N}_{2}\widetilde{\mathcal{B}}^{0}\left({\hat{p}}\right)\\\\[4.73611pt]
-\mathcal{N}^{3}\left(f\right)\\\\[4.73611pt]
\mathcal{B}^{2}\left(\hat{u}\right)\end{array}\right]\;.$ (40)

In (40), we see that the matrices $\mathbb{E}^{3,2}$, $\mathbb{N}_{2}$, are
sparse and metric-free. By metric-free, we mean independent of the size of the
elements, the shape of the elements (orthogonal or highly curved) or the order
of the approximation. All the metric dependence is contained in the basis
functions and therefore in the mass matrix
$\mathbb{M}_{\mathbb{K}^{-1}}^{(2)}$. Using static condensation (40) can be
solved efficiently for the trace variables
$\widetilde{\mathcal{B}}^{0}\left({\lambda}\right)$ only. We first solve for a
global system of $\widetilde{\mathcal{B}}^{0}\left({\lambda}\right)$ by

$\left(\mathbb{N}_{2}^{\intercal}\mathbb{A}\mathbb{N}_{2}\right)\widetilde{\mathcal{B}}^{0}\left({\lambda}\right)=-\mathcal{B}^{2}\left(\hat{u}\right)+\mathbb{N}_{2}^{\intercal}\left({\mathbb{M}_{\mathbb{K}^{-1}}^{(2)}}^{-1}{\mathbb{E}^{3,2}}^{\intercal}\left({\mathbb{E}^{3,2}}{\mathbb{M}_{\mathbb{K}^{-1}}^{(2)}}^{-1}{\mathbb{E}^{3,2}}^{\intercal}\right)^{-1}\mathcal{N}^{3}\left(f\right)-\mathbb{A}\mathbb{N}_{2}\widetilde{\mathcal{B}}^{0}\left({\hat{p}}\right)\right)$
(41)

where,

$\mathbb{A}={\mathbb{M}_{\mathbb{K}^{-1}}^{(2)}}^{-1}-{\mathbb{M}_{\mathbb{K}^{-1}}^{(2)}}^{-1}{\mathbb{E}^{3,2}}^{\intercal}\left({\mathbb{E}^{3,2}}{\mathbb{M}_{\mathbb{K}^{-1}}^{(2)}}^{-1}{\mathbb{E}^{3,2}}^{\intercal}\right)^{-1}{\mathbb{E}^{3,2}}{\mathbb{M}_{\mathbb{K}^{-1}}^{(2)}}^{-1}\;.$

The matrix $\mathbb{A}$ is block diagonal and can be constructed efficiently
by evaluating for each sub domain separately as

$\mathbb{A}_{i}={\mathbb{M}_{\mathbb{K}^{-1},i}^{(2)}}^{-1}-{\mathbb{M}_{\mathbb{K}^{-1},i}^{(2)}}^{-1}{\mathbb{E}_{i}^{3,2}}^{\intercal}\left({\mathbb{E}_{i}^{3,2}}{\mathbb{M}_{\mathbb{K}^{-1},i}^{(2)}}^{-1}{\mathbb{E}_{i}^{3,2}}^{\intercal}\right)^{-1}{\mathbb{E}_{i}^{3,2}}{\mathbb{M}_{\mathbb{K}^{-1},i}^{(2)}}^{-1}\;.$
(42)

The inverse matrices in (42) are symmetric, positive definite and are
evaluated using Cholesky decomposition. The left hand side of (41) is also
symmetric, positive definite and the $\lambda$ system is solved using Cholesky
decomposition. The local degrees of freedom in sub domains $\mathcal{M}_{i}$,
$i=1,\ldots,T$ are then evaluated as

$\displaystyle\widetilde{\mathcal{N}}^{0}_{i}\left({p}\right)$
$\displaystyle=$
$\displaystyle\left(\mathbb{E}^{3,2}_{i}{\mathbb{M}_{\mathbb{K}^{-1},i}^{(2)}}^{-1}{\mathbb{E}^{3,2}_{i}}^{\intercal}\right)^{-1}\left(\mathcal{N}_{i}^{3}\left(f\right)+\mathbb{E}^{3,2}_{i}{\mathbb{M}_{\mathbb{K}^{-1},i}^{(2)}}^{-1}{\mathbb{N}_{2}}_{i}\left({\widetilde{\mathcal{B}}_{i}^{0}\left({\hat{p}}\right)+\widetilde{\mathcal{B}}_{i}^{0}\left({\lambda}\right)}\right)\right)$
(43) $\displaystyle\mathcal{N}^{2}_{i}\left(\bm{u}\right)$ $\displaystyle=$
$\displaystyle{\mathbb{M}_{\mathbb{K}^{-1},i}^{(2)}}^{-1}\left({\mathbb{E}^{3,2}_{i}}^{\intercal}\widetilde{\mathcal{N}}_{i}^{0}\left({p}\right)-{\mathbb{N}_{2}}_{i}\left(\widetilde{\mathcal{B}}^{0}\left(\hat{p}\right)_{i}+\widetilde{\mathcal{B}}_{i}^{0}\left({\lambda}\right)\right)\right)\;.$
(44)

The inverse terms in (43) and (44) are already evaluated in (42), and
therefore evaluation of expansion coefficients of velocity and pressure field
is simply a matrix multiplication step.

##  5 Test cases

In this section we present the computational results for two test cases using
the DD formulation. All the simulations are executed on MATLAB release 2020b
on a Macintosh machine with 2.6 GHz Intel Core i7 processor using a single
processor.

###  5.1 Test case I: Manufactured solution

In this section we will solve a test problem from [14] and compare the results
of DD formulation with the continuous formulation for hexahedral elements of
order $N=1,2,3$ with varying mesh refinements. We will show that i) the
results from both the formulation are same up to machine precision, ii) that
the DD formulation has optimal convergence rates, iii) the speed-up in
simulation run times using DD formulation, and iv) comparison of condition
number of the (35) and (41).

Figure 2: Test domain with total number of elements = $3\times 3\times 3$.
Left: reference domain, $\widehat{\Omega}$, Right: domain of the test case,
$\Omega$.

The domain for the problem is obtained by mapping the reference domain
$\left(\xi\;,\eta\;,\zeta\right)\in\widehat{\Omega}=[-1,1]^{3}$ by

$\left\\{\begin{array}[]{l}x=\hat{x}+0.03\cos\left(3\pi\hat{x}\right)\cos\left(3\pi\hat{y}\right)\cos\left(3\pi\hat{z}\right)\\\\[6.45831pt]
y=\hat{y}-0.04\cos\left(3\pi\hat{x}\right)\cos\left(3\pi\hat{y}\right)\cos\left(3\pi\hat{z}\right)\\\\[6.45831pt]
z=\hat{z}+0.05\cos\left(3\pi\hat{x}\right)\cos\left(3\pi\hat{y}\right)\cos\left(3\pi\hat{z}\right)\end{array}\right.\quad\mbox{where,}\quad\left\\{\begin{array}[]{l}\hat{x}=0.5\left(1+\xi\right)\\\\[6.45831pt]
\hat{y}=0.5\left(1+\eta\right)\\\\[6.45831pt]
\hat{z}=0.5\left(1+\zeta\right)\end{array}\right.\;.$ (45)

In Figure 2, in the left plot we show the reference domain $\widehat{\Omega}$
and on the right plot we show the domain of the problem which is obtained
using (45).

The permeability tensor, $\mathbb{K}$, and the exact solution, $p_{ex}$, are
given by

$\mathbb{K}\left(x,y,z\right)=\left[\begin{array}[]{ccc}x^{2}+y^{2}+1&0&0\\\\[4.73611pt]
0&z^{2}+1&\sin\left(xy\right)\\\\[4.73611pt]
0&\sin\left(xy\right)&x^{2}y^{2}+1\end{array}\right]\;,\qquad\mbox{and}\qquad
p_{ex}\left(x,y,x\right)=x+y+z-1.5\;,$

and the right hand side term is given by

$f_{ex}=-\mathrm{div}\ \mathbb{K}\mathrm{grad}p_{ex}\;.$

As in [14], we impose the Dirchlet boundary conditions $\hat{p}$ at
$\hat{x}=0,1$ faces and Neumann boundary conditions $\hat{u}$ at
$\hat{y}=0,1,\hat{z}=0,1$ faces, where

$\hat{p}=p_{ex}|_{\Gamma_{D}}\qquad\mbox{and}\qquad\hat{u}=\left(-\mathbb{K}\mathrm{grad}p_{ex}\cdot\bm{n}\right)|_{\Gamma_{N}}\;.$

For DD formulation of this test case we decompose the domain $\Omega$ with
equal number of sub domains, $K1$, in each direction. For $N=1,2$ cases each
of the sub domains is discretized into $2\times 2\times 2$ elements. For $N=3$
case each sub domain consists of a single element only. Therefore, total
number of elements in one direction, $K$, are $2K1$ for $N=1,2$ case, and $K1$
for $N=3$ case. We define $\widehat{h}=2/K$ as the size of a non-deformed
element of the reference domain $\widehat{\Omega}=[-1,1]^{3}$.

Figure 3: Comparison of pressure and velocity profiles for the continuous
formulation and the DD formulation for $K=3,N=2$ case. In the first column we
see the results from continuous element formulation, in the second column we
see the results from DD formulation, and in the third column we present the
difference in results between both the formulations. In the first row we see
the pressure, in the second row the $x$-component of the velocity field, in
the third row the $y$-component of the velocity field and in the fourth row
the $z$-component of the velocity field.

In Figure 3 we compare the pressure and velocity profiles of the continuous
formulation and the DD formulation for $K=3,N=2$ case. In the first row we
plot the pressure, in the second row we plot the $x$-component of the
velocity, in the third row we plot the $y$-component of the velocity, and in
the fourth row we plot the $z$-component of the velocity. In the first column
we plot the results from continuous element formulation, in the second column
we plot the results from DD formulation and in the third column we plot the
difference of results between both the formulations. In the third column we
see that the maximum difference between the continuous and the DD formulation
for any of the pressure or velocity profiles is $10^{-13}$.

Figure 4: Error convergence results for the DD formulation. Top left: Error
$(\mathrm{div}\,\bm{u}-f_{ex})$ in $L^{2}\left(\Omega\right)$ norm. Top right:
Error $(\bm{u}-\bm{u}_{ex})$ in $H\left(\mathrm{div};\Omega\right)$ norm.
Bottom centre: Error $(p-p_{ex})$ in $H^{1}\left(\Omega\right)$ norm.

In Figure 4 we show the error convergence results from the DD formulation. At
the top-left we show the convergence of the $L^{2}$-error for the constraint
$(\mathrm{div}\ \bm{u}-f_{ex})$, at top-right we show the convergence of the
error in the $H\left(\mathrm{div};\Omega\right)$ norm for the velocity field,
and in the bottom-centre we show the convergence of the error in the
$H^{1}\left(\Omega\right)$ norm (using (31)), for the pressure field. On the
$x$-axis we have the length of the non-deformed element $\widehat{h}$. All the
error plots show optimal rate of convergence of order
$\mathcal{O}\left(N\right)$.

Table 1: Average computational time (in seconds) for set-up and solution of continuous elements formulation. | $K$ | Set-up time | Solve time | Total time |
---|---|---|---|---|---
$N=1$ |  |  |  |  |
| 4 | $<0.1$ | $<0.1$ | $<0.2$ |
| 8 | $<0.1$ | $<0.1$ | $<0.2$ |
| 16 | 0.3 | 6.5 | 6.8 |
| 32 | 3.3 | 1327.1 | 1330.4 |
| 64 | Out of memory |  |  |
$N=2$ |  |  |  |  |
| 4 | $<0.1$ | $<0.1$ | $<0.2$ |
| 8 | $<0.1$ | 4.9 | 5.0 |
| 16 | 1.3 | 1117.3 | 1118.6 |
| 32 | Out of memory |  |  |
$N=3$ |  |  |  |  |
| 4 | $<0.1$ | 0.8 | 0.9 |
| 8 | 0.5 | 96.1 | 96.6 |
| 16 | Out of memory |  |  |
Table 2: Average computational time (in seconds) for set-up and solution of DD formulation. | $K$ | $K1$ | $K2$ | Set-up time | Solve (41) | Solve (44) & (43) | Total | % to solve (41)
---|---|---|---|---|---|---|---|---
$N=1$ |  |  |  |  |  |  |  |
| 4 | 2 | 2 | 0.1 | $<0.1$ | $<0.1$ | $<0.3$ | -
| 8 | 4 | 2 | 0.1 | $<0.1$ | $<0.1$ | $<0.3$ | -
| 16 | 8 | 2 | 0.8 | $<0.1$ | $<0.1$ | $<1.0$ | -
| 32 | 16 | 2 | 6.4 | 0.4 | 0.1 | 6.9 | 5.8 %
| 64 | 32 | 2 | 50.1 | 21.3 | 0.5 | 71.9 | 29.6 %
| 128 | 64 | 2 | 410.3 | 2149.2 | 19.7 | 2579.2 | 83.3 %
$N=2$ |  |  |  |  |  |  |  |
| 4 | 2 | 2 | 0.2 | $<0.1$ | $<0.1$ | $<0.4$ | -
| 8 | 4 | 2 | 0.2 | $<0.1$ | $<0.1$ | $<0.4$ | -
| 16 | 8 | 2 | 1.9 | 0.4 | 0.1 | 2.3 | 13.0 %
| 32 | 16 | 2 | 15.5 | 11.6 | 0.7 | 27.9 | 41.9 %
| 64 | 32 | 2 | 160.9 | 1763.5 | 63.0 | 1987.4 | 88.7 %
$N=3$ |  |  |  |  |  |  |  |
| 4 | 4 | 1 | 0.2 | $<0.1$ | $<0.1$ | $<0.4$ | -
| 8 | 8 | 1 | 0.7 | 0.1 | $<0.1$ | $<0.9$ | 11.1 %
| 16 | 16 | 1 | 5.3 | 2.7 | 0.2 | 8.2 | 32.9 %
| 32 | 32 | 1 | 51.0 | 267.0 | 9.5 | 327.5 | 81.5 %
Table 3: Average computational time (in seconds) for i) continuous formulation, ii) DD formulation, and the speed-up in simulation time (46). | $K$ | Continuous | DD | Speed-up |
---|---|---|---|---|---
$N=1$ |  |  |  |  |
| 4 | $<0.2$ | $<0.3$ | - |
| 8 | $<0.2$ | $<0.3$ | - |
| 16 | 6.8 | $<1.0$ | 6.8x |
| 32 | 1330.4 | 6.9 | 192.8x |
| 64 | Out of memory | 71.9 | - |
$N=2$ |  |  |  |  |
| 4 | $<0.2$ | $<0.4$ | - |
| 8 | 5.0 | $<0.4$ | 12.5x |
| 16 | 1118.6 | 2.3 | 486.4x |
| 32 | Out of memory | 27.9 | - |
$N=3$ |  |  |  |  |
| 4 | 0.9 | $<0.4$ | - |
| 8 | 96.6 | $<0.9$ | 107.3x |
| 16 | Out of memory | 8.2 | - |

In Table 1, Table 2, Table 3 we give the average simulation times to solve the
continuous formulation, the DD formulation and the comparison between the two
cases respectively. The simulation time (in seconds) is the average for five
simulation runs.

$\|\bm{u}\|_{H\left(div\right);\Omega}^{2}=\|\bm{u}\|^{2}_{L^{2}}+\|\mathrm{div}\bm{u}\|^{2}_{L^{2}}\;.$

$\|\bm{u}-\bm{u}_{ex}\|^{2}=\int_{\Omega}(\bm{u}-\bm{u}_{ex})^{2}\mathrm{d}\Omega$

In Table 1, we give the average simulation time to solve for continuous
formulation for elements of order $N=1,2,3$. In the first column we have the
number of elements in one direction, in the second column we have the time
taken to set-up the matrices and in the third column the time taken to solve
(35) and (36). The simulation times that are less than $0.1$ second, have not
been measured with further accuracy. For $N=1$ the system run out of memory
for $K=64$, for $N=2$ at $K=32$, and for $N=3$ at $K=16$.

In Table 2, we give the average simulation time to solve for DD formulation
for elements of order $N=1,2,3$. In the first column, we give the total numer
of elements in one direction. In the second column we give the number of sub
domains in one direction. In the third column we give the number of elements
in each direction within each sub domain. In the fourth column we give the
time to set-up the matrices. In the fifth column we give the time to solve
(41). In the sixth column we give the time to solve (43) and (44). In the
seventh column we give the total time, and in the eigth column we give the %
amount of time spent to evaluate (41). For the cases where simulation times
are less than 0.1 second the accurate simulation times are not determined and
therefore we do not calculate the last column. We observe in the last column
that for all the three cases i.e. $N=1,2,3$, as we increase the total number
of elements the percentage of time spent on solution of the Lagrange
multiplier system increases. In the last case listed in the table for
$N=1,2,3$, the time spent in solution of (41) is above $80\%$. Given that DD
formulation is for large simulations, we see that in this case the majority of
the time is spent on solution of the Lagrange multiplier system (41).

In Table 3 we compare the total time taken to solve the continuous formulation
and the DD formulation for elements of order $N=1,2,3$. The comparison is made
of cases with same total number of elements in one direction, such that
topology of the mesh discretization remains the same. In the first column we
have the total number of elements in one direction $K$, in the second column
we have the total simulation time from continuous elements formulation, in the
third column we have the total simulation time from DD formulation. In the
fourth column we give the speed-up factor as

$\mbox{speed-up}=\frac{simulation\ time\ continuous\ formulation}{simulation\
time\ DD\ formulation}\;.$ (46)

Firstly, we see that for all $N$, continuous formulation runs out of memory
for lower values of $K$. That is, the DD formulation requires relatively less
memory than the continuous formulation. Secondly, the simulation time
required, for the same refinement, is less for DD method than for continuous
formulation in all cases. Thirdly, the speed-up factor increases as the mesh
is more refined, i.e. the total number of elements are increased.

Figure 5: Comparison of condition number of (35) and (41) for $N=1,2,3$.

In Figure 5 we compare the condition number of the global system of (35) and
(41). We observe that in both the cases, for $N=1,2,3$, and same
discretization, the condition number of the Lagrange multiplier system (41) is
higher than the condition number of the pressure unknowns of (35). Moreover,
the rate of growth of condition number for $N=1,2,3$ is similar in both the
cases.

###  5.2 Test case II: SPE 10

Figure 6: Natural logarithm of permeability field for SPE10 case.

In this section we show the results for solution of SPE10 benchmark problem
[15]. It is often used for validation of numerical schemes for reservoir
modelling applications because of its challenging permeability field. The
domain of the problem and the natural logarithm of the permeability field is
shown in Figure 6. The size of the domain is $1200ft\times 2200ft\times
170ft$. The domain is divided into equal blocks of size $20ft\times 10ft\times
2ft$, i.e. $60\times 220\times 85=1122000$ blocks. Each block has a constant
isotropic permeability tensor $\mathbb{K}_{i}=k_{i}\mathbb{I}$. The top $70ft$
of the domain represents the Tarbert formation that has smooth changes in the
permeability field component, and the bottom $100ft$ represents the Upper Ness
formation that has sharp changes in the permeability field component.

The right hand side term is $f_{ex}=0$. We impose Neumann boundary conditions,
$\hat{u}=\bm{u}\cdot\bm{n}=0$, on $\hat{y},\hat{z}=0,1$, and the Dirichlet
boundary conditions, $\hat{p}=1$ on $\hat{x}=0$, and $\hat{p}=0$ on
$\hat{x}=1$.

The numerical solution of this problem using the continuous formulation was
not possible because the system ran out of memory and therefore we only
present the results with DD method. This shows that DD method can also be
advantageous in cases to reduce the memory requirements for large simulations.

Figure 7: The contour plots of pressure field for SPE10 benchmark problem. Left: results from $Case\ 1$. Centre: results from $Case\ 2$. Right: difference between results from $Case\ 1$ and $Case\ 2$. Table 4: Average computational time (in seconds) for set-up and solution of SPE10 case using DD formulation. | $K$ | $K1$ | $K2$ | Set-up time | Solve (41) | Solve (44) & (43) | Total | % (41)
---|---|---|---|---|---|---|---|---
Case 1 | 1122000 | $15\times 55\times 17$ | $4\times 4\times 5$ | 154.0 | 259.4 | 20.4 | 433.8 | 59.8 %
Case 2 | 1122000 | $12\times 44\times 17$ | $5\times 5\times 5$ | 184.2 | 159.5 | 24.3 | 368.0 | 43.3 %

For numerical solution using the DD method, we use two different approaches.
In $Case\ 1$, we divide the domain into $15\times 55\times 17$ sub domains,
and each sub domain is further divided into $4\times 4\times 5$ elements. In
$Case\ 2$, we divide the domain into $12\times 44\times 17$ sub domains, and
each sub domain is further divided into $5\times 5\times 5$ elements. We also
considered a third case with $6\times 22\times 17$ sub domains, and each sub
domain divided into $10\times 10\times 5$ elements, but this case exceeded the
memory bounds of the system. The total number and topological configuration of
the elements remain the same in all the cases. We use only the lowest order
elements, i.e. $N=1$, in both the cases. The results of the pressure field
from both the cases are shown in Figure 7. In the left plot we see the results
from $Case\ 1$, in the centre plot we see the results from $Case\ 2$, and in
the right plot we see the difference between the two results. We observe that
the maximum difference in results is of order $10^{-12}$.

In Table 4 we give the average run time of five simulation runs to solve the
SPE10 case using DD formulation. In the second column we have the total number
of elements, in the third column we have the total number of sub domains, in
the fourth column we have the total number of elements within each sub domain,
in the fifth column we have the set-up time of matrices, in the sixth column
we have the time taken to solve (41), in the seventh column we have the time
taken to solve $p$ & $\bm{u}$. (44) & (43), in the eighth column we have the
total time, and in the ninth column we have the % of time taken to solve the
Lagrange multiplier system (41) with respect to the total time. The results
show that for the same total number and topology of elements, the choice of
sub domain decomposition also affects the simulation times.

This test case demonstrates that use of algebraic dual representations with DD
method can be applied to practical applications.

##  6 Conclusions

In this paper we have presented the use of algebraic dual spaces for DD method
for Darcy flow. We have defined the broken Sobolev spaces and their finite
dimensional counterparts. We have also defined the global finite dimensional
trace space that connects the broken spaces. These spaces are used to solve
the DD formulation of Darcy equations. It is shown that using algebraic dual
representations the matrix representation of continuity constraint across the
sub domains becomes sparse and metric-free. The first test case is a
manufactured solution where it is shown that i) the results from continuous
formulation and DD formulation are precisely the same and that DD formulation
also has optimal rate of convergence of errors, ii) the DD formulation is more
memory efficient, iii) the DD formulation is more efficient in terms of
simulation run times. In the second test case we show that DD scheme can be
used to reduce the memory requirements and solve for large practical
applications.

We have demonstrated that algebraic dual spaces can be used with DD schemes.
In future this work will be extended to broken $H^{1}\left(\Omega\right)$ and
$H\left(\mathrm{curl};\Omega\right)$ spaces, to address problems such as
vector Laplacian, Stokes flow, and linear elasticity.

##  7 Acknowledgements

The authors would like to thank Dr. Hadi Hajibeygi for the processed
permeability data for SPE10 benchmark case.

## References

  * [1] P. A. Raviart, J. M. Thomas, Primal hybrid finite element methods for 2nd order elliptic equations, Mathematics of Computation 31 (1977) 391–413.
  * [2] B. Wohlmuth, A mortar finite element method using dual spaces for the Lagrange multiplier, SIAM Journal of Numerical Analysis (2000) 989–1012.
  * [3] C. Farhat, J. Mandel, F. X. Roux, Optimal convergence properties of the FETI domain decomposition method, Computer Methods in Applied Mechanics and Engineering 115 (1994) 365–385.
  * [4] C. Farhat, M. Lesoinne, P. LeTallec, K. Pierson, D. Rixen, FETI-DP: a dual–primal unified FETI method—part i: A faster alternative to the two-level FETI method, International Journal for Numerical Methods in Engineering 50 (7) (2001) 1523–1544.
  * [5] L. Demkowicz, J. Gopalakrishnan, A class of discontinuous Petrov–Galerkin methods. Part i: The transport equation, Computer Methods in Applied Mechanics and Engineering 199 (23-24) (2010) 1558–1572.
  * [6] L. Demkowicz, J. Gopalakrishnan, A class of discontinuous Petrov–Galerkin methods. ii. Optimal test functions, Numerical Methods for Partial Differential Equations 27 (1) (2011) 70–105.
  * [7] C. Carstensen, L. Demkowicz, J. Gopalakrishnan, Breaking spaces and forms for the DPG method and applications including Maxwell equations, Computers & Mathematics with Applications 72 (2016) 494–522.
  * [8] V. Agoshkov, Poincaré–Steklov’s operators and domain decomposition methods in finite-dimensional spaces, First International Symposium on Domain Decomposition Methods for Partial Differential Equations, SIAM (1988) 73 – 112\.
  * [9] A. Quarteroni, A. Valli, Theory and application of Steklov-Poincaré operators for boundary-value problems, Mathematics and Its Applications 56, Springer (1991) 179–203.
  * [10] B. Cockburn, Static condensation, hybridization, and the devising of the HDG methods, Lecture Notes in Computational Science and Engineering 114, Springer (2015) 129–177.
  * [11] D. Pietro, A. Ern, S. Lemaire, A review of hybrid high order methods: Formulations, computational aspects, comparison with other methods, Lecture Notes in Computational Science and Engineering 114, Springer (2016) 205 – 236\.
  * [12] J. P. M. Almeida, E. A. Maunder, Equilibrium finite element formulations, John Wiley & Sons Ltd (2016).
  * [13] V. Jain, Y. Zhang, A. Palha, M. Gerritsma, Construction and application of algebraic dual polynomial representations for finite element methods on quadrilateral and hexahedral meshes, Computers & Mathematics with Applications 95 (2021) 101–142.
  * [14] M. F. Wheeler, G. Xue, I. Yotov, A multiscale mortar multipoint flux mixed finite element method, ESAIM: Mathematical Modelling and Numerical Analysis 46 (2012) 759–796.
  * [15] SPE comparative solution project.
URL https://www.spe.org/web/csp/datasets/set02.htm