# A new numerical method for div-curl Systems with Low Regularity Assumptions

Shuhao Cao  Department of Mathematics and Statistics, Washington University in
St. Louis, St. Louis, MO 63130<EMAIL_ADDRESS>,  Chunmei Wang  Department of
Mathematics & Statistics, Texas Tech University, Lubbock, TX 79409, USA
<EMAIL_ADDRESS>and  Junping Wang  Division of Mathematical Sciences,
National Science Foundation, Alexandria, VA 22314<EMAIL_ADDRESS>
###### Abstract.

This paper presents a numerical method for div-curl systems with normal
boundary conditions by using a finite element technique known as primal-dual
weak Galerkin (PDWG). The PDWG finite element scheme for the div-curl system
has two prominent features in that it offers not only an accurate and reliable
numerical solution to the div-curl system under the low
$H^{\alpha}$-regularity ($\alpha>0$) assumption for the true solution, but
also an effective approximation of normal harmonic vector fields regardless
the topology of the domain. Results of seven numerical experiments are
presented to demonstrate the performance of the PDWG algorithm, including one
example on the computation of discrete normal harmonic vector fields.

###### Key words and phrases:

finite element methods, weak Galerkin methods, primal-dual weak Galerkin, div-
curl system

###### 2020 Mathematics Subject Classification:

Primary 65N30, 35Q60, 65N12; Secondary 35F45, 35Q61

The research of Shuhao Cao was partially supported by National Science
Foundation Award DMS-1913080.

The research of Chunmei Wang was partially supported by National Science
Foundation Award DMS-1849483.

The research of Junping Wang was supported by the NSF IR/D program, while
working at National Science Foundation. However, any opinion, finding, and
conclusions or recommendations expressed in this material are those of the
author and do not necessarily reflect the views of the National Science
Foundation.

##  1\. Introduction

In this paper we are concerned with the development of new numerical methods
for div-curl systems equipped with normal boundary conditions. For simplicity,
consider the model problem that seeks a vector field $\bm{u}=\bm{u}({\bm{x}})$
satisfying

(1.1) $\displaystyle\nabla\cdot(\varepsilon\bm{u})$
$\displaystyle=f,\qquad{\rm in}\ \Omega,$ (1.2)
$\displaystyle\nabla\times\bm{u}$ $\displaystyle=\bm{g},\qquad{\rm in}\
\Omega,$ (1.3) $\displaystyle\varepsilon\bm{u}\cdot\bm{n}$
$\displaystyle=\phi_{1},\qquad{\rm on}\ \Gamma,$

where $\Omega\subset{\mathbb{R}}^{3}$ is a polyhedral domain that is bounded,
open, and connected. $\Gamma=\partial\Omega$ is the boundary of $\Omega$ and
is assumed to be the union of a finite number of disjoint surfaces
$\Gamma=\bigcup_{i=0}^{L}\Gamma_{i}$, where $\Gamma_{0}$ is the exterior
boundary of $\Omega$, $\Gamma_{i}\ (i=1,\cdots,L)$ are the other connected
components with finite surface areas, and $L$ denotes the number of holes in
the domain $\Omega$ geometrically and is known as the second Betti number of
$\Omega$ or the dimension of the second de Rham cohomology group of $\Omega$.
The Lebesgue-integrable real-valued function $f=f({\bm{x}})$ and the vector
field $\bm{g}=\bm{g}({\bm{x}})$ are given in the domain $\Omega$, the
coefficient matrix $\varepsilon=\\{\varepsilon_{ij}({\bm{x}})\\}_{3\times 3}$
is symmetric and uniformly positive definite in $\Omega$, and the entries
$\varepsilon_{ij}({\bm{x}})$ ($i,j=1,\cdots,3$) are in $L^{\infty}(\Omega)$.
The normal boundary data $\phi_{1}$ is a given distribution in
$H^{-\frac{1}{2}}(\Gamma)$.

The solution uniqueness for the normal boundary value problem (1.1)-(1.3)
depends on the topology of the domain $\Omega$. It is well-known that the
solution uniqueness holds true for simply connected $\Omega$, while the
solution is unique up to a normal $\varepsilon$-harmonic function in
$\mathbb{H}_{\varepsilon n,0}(\Omega)$ defined in (2.2) when the domain
$\Omega$ is not simply connected. The dimension of $\mathbb{H}_{\varepsilon
n,0}(\Omega)$ is identical to the first Betti number of $\Omega$, which is the
rank of the first homology group of $\Omega$ or the dimension of the first de
Rham cohomology group of $\Omega$.

The div-curl system (1.1)-(1.2) arises in many applications such as
electromagnetic fields and fluid mechanics. Computational electro-magnetics
plays an important role in many areas of science and engineering such as
radar, satellite, antenna design, waveguides, optical fibers, medical imaging
and design of invisible cloaking devices [14]. In linear magnetic fields, the
function $f({\bm{x}})$ vanishes, ${\bm{u}}$ represents the magnetic field
intensity and $\varepsilon({\bm{x}})$ is the inverse of the magnetic
permeability tensor. In fluid mechanics fields, the coefficient matrix
$\varepsilon({\bm{x}})$ is diagonal with diagonal entries being the local mass
density. In electrostatics fields, $\varepsilon({\bm{x}})$ is the permittivity
matrix.

Several numerical methods based on finite element approaches have been
proposed and analyzed for the div-curl system (1.1)-(1.2). For example, a
covolume method was developed in three dimensional space by using the Voronoi-
Delaunay mesh pairs [18]. In [5], a finite element method based on different
test and trial spaces was analyzed for a div-curl system. In [16], the authors
introduced a mimetic finite difference scheme for the three dimensional
magneto-static problems on general polyhedral partitions. In [10], the authors
developed a mixed finite element method for three dimensional axisymmetric
div-curl systems through a dimension reduction technique based on the
cylindrical coordinates in simply connected and axisymmetric domains. In [4],
the authors developed a least-squares finite element method for two types
boundary value problems. Another least-squares method was proposed for the
div-curl problem based on discontinuous elements on nonconvex polyhedral
domains in [3]. In [21], the authors proposed a weak Galerkin finite element
method for the div-curl system with either normal or tangential boundary
conditions. Another weak Galerkin scheme was introduced in [15] by using a
least-squares approach for the div-curl problem. In [17], the authors
developed a primal-dual weak Galerkin finite element method for the div-curl
system with tangential boundary conditions and proved that the scheme works
well with low-regularity assumptions on the exact solution.

Two major challenges for the div-curl model problem (1.1)-(1.3) are the low-
regularity of the true solution ${\bm{u}}$ and the difficulties in
approximating the normal $\varepsilon$-harmonic vector space
$\mathbb{H}_{\varepsilon n,0}(\Omega)$ due to the topological complexity of
the domain $\Omega$. The goal of this paper is to address both challenges
through a new numerical method based on the primal-dual weak Galerkin (PDWG)
finite element aproach originated in [23] and further developed in [20, 25,
26, 7, 22, 13] for various model problems. Our PDWG numerical method for
(1.1)-(1.3) has two prominent features over the existing numerical methods:
(1) it offers an effective approximation of the normal $\varepsilon$-harmonic
vector space $\mathbb{H}_{\varepsilon n,0}(\Omega)$ regardless of the topology
of the domain $\Omega$; and (2) it provides an accurate and reliable numerical
solution for the div-curl system (1.1)-(1.3) with low $H^{\alpha}$-regularity
($\alpha>0$) assumption for the true solution ${\bm{u}}$.

The paper is organized as follows. Section 2 is devoted to notations and the
derivation of a weak formulation for the div-curl system (1.1)-(1.3) that
involves no partial derivatives over the vector field ${\bm{u}}$. Section 3
offers a brief review on the discrete weak gradient and discrete weak curl
operators. Section 4 is dedicated to the presentation of the PDWG algorithm
for the div-curl problem, together with an algorithm for computing discrete
normal $\varepsilon$-harmonic vector fields. Section 5 is devoted to a
discussion of the solution existence and uniqueness for the PDWG scheme.
Section 6 contains a convergence theory for the PDWG approximation, and
Section 7 demonstrates the performance of the PDWG algorithm through seven
test examples.

##  2\. Notations and Preliminaries

We follow the usual notation for Sobolev spaces and norms, see for example [9,
12]. For an open bounded domain $D\subset{\mathbb{R}}^{3}$ with Lipschitz
continuous boundary and any given real number $s\geq 0$, we use
$\|\cdot\|_{s,D}$ and $|\cdot|_{s,D}$ to denote the norm and seminorm in the
Sobolev space $H^{s}(D)$, respectively. The space $H^{0}(D)$ coincides with
$L^{2}(D)$, for which the norm and the inner product are denoted by
$\|\cdot\|_{D}$ and $(\cdot,\cdot)_{D}$, respectively. We use
$H(div_{\varepsilon};D)$ to denote the closed subspace of $[L^{2}(D)]^{2}$ so
that $\nabla\cdot(\varepsilon{\bm{v}})\in L^{2}(D)$. The space $H(div;D)$
corresponds to the case of $\varepsilon=I$. Analogously, we use $H(curl;D)$ to
denote the closed subspace of $[L^{2}(D)]^{2}$ so that
$\nabla\times{\bm{v}}\in[L^{2}(D)]^{3}$. Denote by

$H_{0}(curl;D):=\\{{\bm{v}}\in H(curl;D),\ {\bm{v}}\times{\bm{n}}=0\mbox{ on
}\partial D\\}$

the closed subspace with vanishing tangential boundary values. When
$D=\Omega$, we shall drop the script $D$ in the notations. For simplicity, we
shall use “$\lesssim$ ” to denote “less than or equal to up to a general
constant independent of the mesh size or functions appearing in the
inequality”.

Introduce the following Sobolev space

(2.1) ${\mathbb{W}_{\varepsilon}(\Omega)}=\\{{\bm{v}}\in H_{0}(curl)\cap
H(div_{\varepsilon}),\ \nabla\cdot(\varepsilon{\bm{v}})=0,\
\langle\varepsilon{\bm{v}}\cdot{\bm{n}}_{i},1\rangle_{\Gamma_{i}}=0,\
i=1,\cdots,L\\}.$

A vector field ${\bm{v}}\in[L^{2}(\Omega)]^{3}$ is said to be
$\varepsilon$-harmonic on $\Omega$ if it is $\varepsilon$-solenoidal and
irrotational on $\Omega$. The space of normal $\varepsilon$-harmonic vector
fields, denoted by $\mathbb{H}_{\varepsilon n,0}(\Omega)$, consists of all
$\varepsilon$-harmonic vector fields satisfying the zero normal boundary
condition; i.e.,

(2.2) $\mathbb{H}_{\varepsilon
n,0}(\Omega)=\\{{\bm{v}}\in[L^{2}(\Omega)]^{3}:\ {\nabla\times}{\bm{v}}=0,\
\nabla\cdot(\varepsilon{\bm{v}})=0,\ \varepsilon{\bm{v}}\cdot{\bm{n}}=0\mbox{
on }\Gamma\\}.$

When $\varepsilon=I$ is the identity matrix, the space
$\mathbb{H}_{\varepsilon n,0}(\Omega)$ shall be denoted as
$\mathbb{H}_{n,0}(\Omega)$. Analogously, the space of tangential
$\varepsilon$-harmonic vector fields, denoted by
$\mathbb{H}_{\varepsilon\tau,0}(\Omega)$, consists of all
$\varepsilon$-harmonic vector fields satisfying the zero tangential boundary
condition; i.e.,

$\mathbb{H}_{\varepsilon\tau,0}(\Omega)=\\{{\bm{v}}\in[L^{2}(\Omega)]^{3}:\
{\nabla\times}{\bm{v}}=0,\ \nabla\cdot(\varepsilon{\bm{v}})=0,\
{\bm{v}}\times{\bm{n}}=0\mbox{ on }\Gamma\\}.$

###  2.1. A weak formulation

By testing the equation (1.1) against any $\varphi\in H^{1}(\Omega)$ and then
using the normal boundary condition (1.3) we obtain

(2.3)
$(\bm{u},\varepsilon\nabla\varphi)=\langle\phi_{1},\varphi\rangle-(f,\varphi),\qquad\forall\varphi\in
H^{1}(\Omega).$

Next, we test the equation (1.2) against any $\bm{w}\in H_{0}(curl;\Omega)$ to
obtain

(2.4)
$\displaystyle(\bm{u},\nabla\times\bm{w})=(\bm{g},\bm{w}),\qquad\forall\bm{w}\in
H_{0}(curl;\Omega).$

Summing the equations (2.3) and (2.4) gives the following

$\displaystyle(\bm{u},\varepsilon\nabla\varphi+\nabla\times\bm{w})=(\bm{g},\bm{w})-(f,\varphi)+\langle\phi_{1},\varphi\rangle$

for all $\varphi\in H^{1}(\Omega)$ and $\bm{w}\in H_{0}(curl;\Omega)$.

######  Definition 2.1.

A vector-valued function $\bm{u}\in[L^{2}(\Omega)]^{3}$ is said to be a weak
solution of the normal boundary value problem for the div-curl system
(1.1)-(1.3) if it satisfies the following equation

(2.5)
$\displaystyle(\bm{u},\varepsilon\nabla\varphi+\nabla\times{\bm{\psi}})=(\bm{g},{\bm{\psi}})-(f,\varphi)+\langle\phi_{1},\varphi\rangle$

for all $\varphi\in H^{1}(\Omega)$ and ${\bm{\psi}}\in H_{0}(curl;\Omega)$.

The solution to the variational problem (2.5) is generally non-unique. In
fact, the homogeneous version of (2.5) seeks ${\bm{u}}\in[L^{2}(\Omega)]^{3}$
satisfying

(2.6)
$\displaystyle(\bm{u},\varepsilon\nabla\varphi+\nabla\times{\bm{\psi}})=0\qquad\forall\varphi\in
H^{1}(\Omega),\ {\bm{\psi}}\in H_{0}(curl;\Omega).$

The equation (2.6) is easily satisfied by any $\varepsilon$-harmonic function
${\bm{u}}=\bm{\eta}\in\mathbb{H}_{\varepsilon n,0}(\Omega)$, and hence the
solution non-uniqueness when the $\varepsilon$-harmonic space
$\mathbb{H}_{\varepsilon n,0}(\Omega)$ has a positive dimension. The solution
to the div-curl system with normal boundary condition is unique when the
solution is further required to be $\varepsilon$-weighted $L^{2}$ orthogonal
to $\mathbb{H}_{\varepsilon n,0}(\Omega)$.

###  2.2. An extended weak formulation

Denote by $H_{0c}^{1}(\Omega)$ the space of functions in $H^{1}(\Omega)$ with
vanishing value on $\Gamma_{0}$ and constant values on other connected
components of the boundary; i.e.,

$H_{0c}^{1}(\Omega)=\\{\phi\in H^{1}(\Omega):\ \phi|_{\Gamma_{0}}=0,\
\phi|_{\Gamma_{i}}=\alpha_{i},\ i=1,\dotsc,L\\}.$

Introduce the following bilinear form:

(2.7)
$B({\bm{u}},s;\varphi,{\bm{\psi}}):=({\bm{u}},\varepsilon\nabla\varphi+\nabla\times{\bm{\psi}})+({\bm{\psi}},\varepsilon\nabla
s).$

The extended weak formulation for the normal boundary value problem of the
div-curl system seeks $({\bm{u}},s)\in[L^{2}(\Omega)]^{3}\times
H_{0c}^{1}(\Omega)$ such that

(2.8)
$B({\bm{u}},s;\varphi,{\bm{\psi}})=F(\varphi,{\bm{\psi}})\qquad\forall\varphi\in
H^{1}(\Omega),\forall{\bm{\psi}}\in H_{0}(curl;\Omega),$

where

$F(\varphi,{\bm{\psi}})=(\bm{g},{\bm{\psi}})-(f,\varphi)+\langle\phi_{1},\varphi\rangle.$

Note that by testing the curl equation in the div-curl system against any
$\rho\in H^{1}_{0c}(\Omega)$ we have

(2.9) $(\bm{g},\nabla\rho)=0,\qquad\forall\rho\in H_{0c}^{1}(\Omega),$

which gives rise to the following compatibility condition:

(2.10)
$\nabla\cdot\bm{g}=0,\quad\langle\bm{g}\cdot{\bm{n}}_{i},1\rangle_{\Gamma_{i}}=0,\
\text{for}\ i=1,\cdots,L.$

######  Theorem 2.1.

Under the compatibility condition (2.10) for $\bm{g}$, the solution
$({\bm{u}};s)$ of (2.8) satisfies the following equations:

(2.11) $\displaystyle\nabla\cdot(\varepsilon{\bm{u}})$ $\displaystyle=$
$\displaystyle f,\qquad\mbox{ in }\Omega,$ (2.12)
$\displaystyle\nabla\times{\bm{u}}$ $\displaystyle=$
$\displaystyle\bm{g},\qquad\mbox{ in }\Omega,$ (2.13) $\displaystyle s$
$\displaystyle=$ $\displaystyle 0,\qquad\mbox{ in }\Omega,$ (2.14)
$\displaystyle\varepsilon{\bm{u}}\cdot{\bm{n}}$ $\displaystyle=$
$\displaystyle\phi_{1},\qquad\mbox{on }\Gamma.$

###### Proof.

By letting $\bm{\psi}=0$ in (2.8) we have

$(\bm{u},\varepsilon\nabla\varphi)=-(f,\varphi)+\langle\phi_{1},\varphi\rangle_{\Gamma}$

for all $\varphi\in H^{1}(\Omega)$. It follows that
$\nabla\cdot(\varepsilon{\bm{u}})=f$ and
$\varepsilon{\bm{u}}\cdot{\bm{n}}=\phi_{1}\mbox{ on }\Gamma$, which leads to
(2.11) and (2.14). Next, by letting $\varphi=0$ in (2.8) we arrive at

$\displaystyle(\bm{u},\nabla\times\bm{\psi})+(\bm{\psi},\varepsilon\nabla
s)=(\bm{g},\bm{\psi}),\qquad\forall\bm{\psi}\in H_{0}(curl;\Omega),$

which leads to

$\displaystyle(\nabla\times\bm{u}+\varepsilon\nabla
s,\bm{\psi})=(\bm{g},\bm{\psi}),$

and thus

(2.15) $\displaystyle\nabla\times\bm{u}+\varepsilon\nabla
s=\bm{g},\qquad\mbox{in}\ \Omega.$

Now from (2.15) we have

$(\nabla\times\bm{u}+\varepsilon\nabla s,\nabla s)=(\bm{g},\nabla s),$

which, by the usual integration by parts, gives

$\langle{\bm{n}}\times\bm{u},\nabla s\rangle_{\Gamma}+(\varepsilon\nabla
s,\nabla s)=(\bm{g},\nabla s),$

and by the boundary condition of $s=const$ on each $\Gamma_{i}$ and the
compatibility condition (2.10)

$(\varepsilon\nabla s,\nabla s)=\langle\bm{u},{\bm{n}}\times\nabla
s\rangle_{\Gamma}+(\bm{g},\nabla s)=0.$

It follows that $\nabla s=0$ so that $s\equiv 0$. This completes the proof of
the theorem. ∎

The homogeneous dual problem of (2.8) seeks $\lambda\in
H^{1}(\Omega)/\mathbb{R}$ and ${\bm{q}}\in H_{0}(curl;\Omega)$ satisfying

(2.16)
$B({\bm{v}},r;\lambda,{\bm{q}})=0\qquad\forall{\bm{v}}\in[L^{2}(\Omega)]^{3},\
r\in H_{0c}^{1}(\Omega).$

######  Theorem 2.2.

The solution to the homogeneous dual problem (2.16) is unique.

###### Proof.

The problem (2.16) can be rewritten as

(2.17)
$(\bm{v},\varepsilon\nabla\lambda+\nabla\times{\bm{q}})+({\bm{q}},\varepsilon\nabla
r)=0$

for all ${\bm{v}}\in[L^{2}(\Omega)]^{3}$ and $r\in H_{0c}^{1}(\Omega)$. Note
that the test against $r\in H^{1}_{0c}(\Omega)$ and $\bm{v}=0$ ensures
$\nabla\cdot(\varepsilon{\bm{q}})=0$ and
$\langle\varepsilon{\bm{q}}\cdot{\bm{n}},1\rangle_{\Gamma_{i}}=0$ for all
$i\in\\{1,\cdots,L\\}$. In addition, by letting $r=0$ and varying
$\bm{v}\in[L^{2}(\Omega)]^{3}$ we arrive at

$\varepsilon\nabla\lambda+\nabla\times{\bm{q}}=0,$

which, by testing against $\nabla\lambda$, leads to

$(\varepsilon\nabla\lambda,\nabla\lambda)=0,$

so that $\lambda\equiv 0$ and hence

$\nabla\times\bm{q}=0.$

Thus, we have

$\bm{q}\in\mathbb{H}_{\varepsilon\tau,0}(\Omega),\quad\langle\varepsilon\bm{q}\cdot{\bm{n}}_{i},1\rangle_{\Gamma_{i}}=0\
\text{for}\ i=1,\cdots,L,$

which yields $\bm{q}\equiv 0$. ∎

##  3\. Discrete Weak Gradient and Weak Curl

The extended weak formulation (2.8) is based on the gradient and curl
differential operators. In this section, we shall review the notion of
discrete weak gradient and weak curl which forms a corner stone of the weak
Galerkin finite element method. To this end, let $T$ be a polyhedral domain
with boundary $\partial T$ and unit outward normal direction $\bm{n}$ on
$\partial T$. Define the space of weak functions in $T$ by

$W(T)=\\{v=\\{v_{0},v_{b}\\}:v_{0}\in L_{2}(T),v_{b}\in L_{2}(\partial T)\\},$

where $v_{0}$ represents the value of $v$ in the interior of $T$, and $v_{b}$
represents some specific boundary information of $v$. Analogously, define the
space of vector-valued weak functions on $T$ by

$V(T)=\\{\bm{v}=\\{\bm{v}_{0},\bm{v}_{b}\\}:\bm{v}_{0}\in[L_{2}(T)]^{3},\bm{v}_{b}\in[L_{2}(\partial
T)]^{3}\\}.$

Let $P_{j}(T)$ the space of polynomials on $T$ with total degree $j$ and less.
For any $v\in W(T)$, the discrete weak gradient, denoted by $\nabla_{w,j,T}v$,
is defined as the unique vector-valued polynomial in $[P_{j}(T)]^{3}$
satisfying

(3.1)
$(\nabla_{w,j,T}v,\bm{\varphi})_{T}=-(v_{0},\nabla\cdot\bm{\varphi})_{T}+\langle
v_{b},\bm{\varphi}\cdot{\bm{n}}\rangle_{\partial
T},\quad\forall\;\bm{\varphi}\in[P_{j}(T)]^{3}.$

Analogously, the discrete weak curl of ${\bm{v}}\in V(T)$, denoted by
$\nabla_{w,j,T}\times{\bm{v}}$, is defined as the unique vector-valued
polynomial in $[P_{j}(T)]^{3}$, satisfying

(3.2)
$(\nabla_{w,j,T}\times{\bm{v}},{\bm{\varphi}})_{T}=({\bm{v}}_{0},\nabla\times{\bm{\varphi}})_{T}-\langle{\bm{v}}_{b}\times{\bm{n}},{\bm{\varphi}}\rangle_{\partial
T},\quad\forall\;\bm{\varphi}\in[P_{j}(T)]^{3}.$

##  4\. PDWG Numerical Algorithm

Let ${\mathcal{T}}_{h}$ be a finite element partition of the domain $\Omega$
consisting of polyhedra that are shape-regular [27]. Denote by
${\mathcal{E}}_{h}$ the set of faces in ${\mathcal{T}}_{h}$ and
${\mathcal{E}}_{h}^{0}={\mathcal{E}}_{h}\setminus\partial\Omega$ the set of
interior faces. Denote by $h_{T}$ the diameter of the element
$T\in{\mathcal{T}}_{h}$ and $h=\max_{T\in{\mathcal{T}}_{h}}h_{T}$ the meshsize
of the partition ${\mathcal{T}}_{h}$.

For a given integer $k\geq 0$, we introduce the following finite element
spaces subordinated to ${\mathcal{T}}_{h}$:

$\begin{split}{\bm{V}}_{h}=&\\{\bm{v}:\ \bm{v}|_{T}\in[P_{k}(T)]^{3},\forall
T\in{\mathcal{T}}_{h}\\},\\\ S_{h}=&\\{\\{s_{0},s_{b}\\}:\ s_{0}|_{T}\in
P_{k}(T),s_{b}|_{\partial T}\in P_{k}({\partial T}),\forall
T\in{\mathcal{T}}_{h},s_{b}|_{\Gamma_{0}}=0,s_{b}|_{\Gamma_{i}}=c_{i}\\},\\\
M_{h}=&\\{\\{\varphi_{0},\varphi_{b}\\}:\ \varphi_{0}|_{T}\in
P_{k}(T),\varphi_{b}|_{\partial T}\in P_{k}({\partial T}),\forall
T\in{\mathcal{T}}_{h},\ \int_{\Omega}\varphi_{0}=0\\},\\\
{\bm{W}}_{h}=&\\{\bm{\psi}=\\{\bm{\psi}_{0},\bm{\psi}_{b}\\}:\
\bm{\psi}_{0}|_{T}\in[P_{k}(T)]^{3},\bm{\psi}_{b}|_{\partial T}\in
G_{k}({\partial T}),\forall
T\in{\mathcal{T}}_{h},\bm{\psi}_{b}|_{\Gamma}=0\\},\end{split}$

where $G_{k}({\partial T}):=[P_{k}(\sigma)]^{3}\times{\bm{n}}_{\sigma}$ is the
space of polynomials of degree $k$ in the tangent space of ${\partial T}$.

For simplicity of notation, for $\sigma\in S_{h}$ or $\sigma\in M_{h}$, denote
by $\nabla_{w}\sigma$ the discrete weak gradient $\nabla_{w,k,T}\sigma$
computed by using (3.1) on each element $T$; i.e.,

$(\nabla_{w}\sigma)|_{T}=\nabla_{w,k,T}(\sigma|_{T}),\qquad\sigma\in S_{h}\
\text{or}\ \sigma\in M_{h}.$

Analogously, for ${\bm{q}}\in{\bm{W}}_{h}$, denote by
$\nabla_{w}\times{\bm{q}}$ the discrete weak curl
$\nabla_{w,k,T}\times{\bm{q}}$ computed by using (3.2) on each element $T$;
i.e.,

$(\nabla_{w}\times{\bm{q}})|_{T}=\nabla_{w,k,T}\times({\bm{q}}|_{T}),\qquad{\bm{q}}\in{\bm{W}}_{h}.$

An approximation of the bilinear form $B(\cdot;\cdot)$ is given as follows:

(4.1)
$B_{h}({\bm{u}}_{h},s_{h};\varphi,{\bm{\psi}}):=({\bm{u}}_{h},\varepsilon\nabla_{w}\varphi+\nabla_{w}\times{\bm{\psi}})+({\bm{\psi}}_{0},\varepsilon\nabla_{w}s_{h})$

for ${\bm{u}}_{h}\in{\bm{V}}_{h},\ s_{h}\in S_{h},\ \varphi\in M_{h},\
{\bm{\psi}}\in{\bm{W}}_{h}$.

The following is the PDWG finite element method for the div-curl model system
(1.1)-(1.3).

######  Algorithm 1 (PDWG Algorithm).

For an approximate solution of (1.1)-(1.3), one may compute
${\bm{u}}_{h}\in{\bm{V}}_{h}$, together with three auxiliary variables
$s_{h}\in S_{h}$, $\lambda_{h}\in M_{h}$, and ${\bm{q}}_{h}\in{\bm{W}}_{h}$
satisfying

(4.2)
$\left\\{\begin{array}[]{rl}s_{1}(\lambda_{h},{\bm{q}}_{h};\varphi,{\bm{\psi}})+B_{h}({\bm{u}}_{h},s_{h};\varphi,{\bm{\psi}})&=F(\varphi,{\bm{\psi}}),\quad\forall\varphi\in
M_{h},\ {\bm{\psi}}\in{\bm{W}}_{h},\\\
-s_{2}(s_{h},r)+B_{h}(\bm{v},r;\lambda_{h},{\bm{q}}_{h})&=0,\quad\forall\bm{v}\in{\bm{V}}_{h},\
r\in S_{h}.\end{array}\right.$

Here the stabilizer $s_{1}$ is given by

$\displaystyle s_{1}(\lambda_{h},{\bm{q}}_{h};\varphi,\bm{\psi})$
$\displaystyle=$
$\displaystyle\rho_{1}\sum_{T\in{\mathcal{T}}_{h}}h_{T}^{-1}\langle\lambda_{0}-\lambda_{b},\varphi_{0}-\varphi_{b}\rangle_{\partial
T}$
$\displaystyle+\rho_{2}\sum_{T\in{\mathcal{T}}_{h}}h_{T}^{-1}\langle\bm{q}_{0}\times\bm{n}-\bm{q}_{b}\times\bm{n},\bm{\psi}_{0}\times\bm{n}-\bm{\psi}_{b}\times\bm{n}\rangle_{\partial
T},$

and $s_{2}$ is defined accordingly in the space $M_{h}$ as follows

$\displaystyle
s_{2}(s_{h};r)=\rho_{3}\sum_{T\in{\mathcal{T}}_{h}}h_{T}^{-\gamma}\langle
s_{0}-s_{b},r_{0}-r_{b}\rangle_{\partial T},$

where $\gamma\geq-1$ and $\rho_{i}>0$ are parameters with values at user’s
discretion.

The PDWG scheme (4.2) further provides an approximation of the space of normal
$\varepsilon$-harmonic vector fields $\mathbb{H}_{\varepsilon n,0}(\Omega)$,
as revealed by Theorem 6.2 in that the difference
$\bm{\eta}_{h}={\mathcal{Q}}_{h}{\bm{u}}-{\bm{u}}_{h}$ is sufficiently close
to a true normal $\varepsilon$-harmonic vector field $\bm{\eta}$. For this
purpose, we introduce the following notation of discrete normal
$\varepsilon$-harmonic functions.

######  Definition 4.1 (discrete normal $\varepsilon$-harmonic functions).

A vector field $\bm{\eta}_{h}\in{\bm{V}}_{h}$ is said to be a discrete normal
$\varepsilon$-harmonic function if there exists a vector field ${\bm{u}}\in
H(div_{\varepsilon};\Omega)\cap H(curl;\Omega)$ such that

(4.4) $\bm{\eta}_{h}={\mathcal{Q}}_{h}{\bm{u}}-{\bm{u}}_{h},$

where ${\mathcal{Q}}_{h}$ is the $L^{2}$ projection operator onto the finite
element space ${\bm{V}}_{h}$ and ${\bm{u}}_{h}$ is the solution of (4.2) for a
div-curl system (1.1)-(1.3) with load functions $f$, $\bm{g}$, and $\phi_{1}$
determined by ${\bm{u}}$.

In practical computation, a discrete normal $\varepsilon$-harmonic function
can be readily obtained from (4.4) by choosing a smooth vector field
${\bm{u}}$ and one solving of the PDWG system (4.2).

##  5\. Solution Existence and Uniqueness

Introduce two semi-norms as follows:

(5.1)
${|||}(\lambda_{h},{\bm{q}}_{h}){|||}=\Big{(}s_{1}(\lambda_{h},{\bm{q}}_{h};\lambda_{h},{\bm{q}}_{h})\Big{)}^{\frac{1}{2}},\qquad\lambda_{h}\in
M_{h},\ {\bm{q}}_{h}\in{\bm{W}}_{h},$ (5.2)
${|||}s_{h}{|||}=\Big{(}s_{2}(s_{h};s_{h})\Big{)}^{\frac{1}{2}},\qquad
s_{h}\in S_{h}.$

For simplicity, assume that $\varepsilon$ is piecewise constant with respect
to the partition ${\mathcal{T}}_{h}$. Note that all the results can be
generalized to piecewise smooth $\varepsilon$ without any difficulty.

Denote by $Q_{0}$ the $L^{2}$ projection operator onto $P_{k}(T)$ and $Q_{b}$
the $L^{2}$ projection operator onto $P_{k}(\sigma)$ on each face
$\sigma\in\partial T$. Denote by $Q_{h}$ the projection operator onto the weak
finite element space $S_{h}$ or $M_{h}$ such that

$(Q_{h}w)|_{T}=\\{Q_{0}w|_{T},Q_{b}w|_{{\partial T}}\\}.$

Analogously, we use $\mathbb{Q}_{0}$, $\mathbb{Q}_{b}$ and $\mathbb{Q}_{h}$ to
denote the $L^{2}$ projection operators onto $[P_{k}(T)]^{3}$,
$G_{k}(\sigma)=[P_{k}(\sigma)]^{3}\times{\bm{n}}_{\sigma}$, and
${\bm{W}}_{h}$, respectively. The $L^{2}$ projection operator onto the finite
element space ${\bm{V}}_{h}$ is denoted as ${\mathcal{Q}}_{h}$.

######  Lemma 5.1.

[27] The $L^{2}$ projections $Q_{h}$ and ${\mathcal{Q}}_{h}$ satisfy the
commutative property

(5.3) $\nabla_{w}(Q_{h}w)={{\mathcal{Q}}}_{h}(\nabla w),\qquad\forall w\in
H^{1}(T),$ (5.4)
$\nabla_{w}\times(\mathbb{Q}_{h}{\bm{\psi}})={\mathcal{Q}}_{h}(\nabla\times{\bm{\psi}}),\qquad\forall{\bm{\psi}}\in
H(curl;T).$

######  Theorem 5.2.

The kernel of the matrix of the PDWG finite element scheme (4.2) is given by

$K_{h}=\\{({\bm{v}}_{h},s_{h}=0,\lambda_{h}=0,{\bm{q}}_{h}=0):\
{\bm{v}}_{h}\in{\bm{V}}_{h}\cap\mathbb{H}_{\varepsilon n,0}(\Omega)\\}.$

In other words, the kernel of the matrix of the PDWG scheme (4.2) is
isomorphic to the subspace of $\mathbb{H}_{\varepsilon n,0}(\Omega)$
consisting of harmonic functions that are piecewise polynomial of degree $k$.

###### Proof.

Let $({\bm{u}}_{h},s_{h},\lambda_{h},\bm{q}_{h})$ be a solution of (4.2) with
homogeneous data. It follows that

(5.5) $\displaystyle
s_{1}(\lambda_{h},\bm{q}_{h};\lambda_{h},\bm{q}_{h})=0,\quad
s_{2}(s_{h},s_{h})=0,$ (5.6)
$\displaystyle(\bm{u}_{h},\varepsilon\nabla_{w}\varphi+\nabla_{w}\times\bm{\psi})+(\bm{\psi}_{0},\varepsilon\nabla_{w}s_{h})=0,\qquad\forall\varphi\in
M_{h},\bm{\psi}\in{\bm{W}}_{h},$ (5.7)
$\displaystyle(\bm{q}_{0},\varepsilon\nabla_{w}r)+(\bm{v},\varepsilon\nabla_{w}\lambda_{h}+\nabla_{w}\times\bm{q}_{h})=0,\qquad\forall\bm{v}\in{\bm{V}}_{h},r\in
S_{h}.$

From (5.5) we have

(5.8)
$\displaystyle\lambda_{0}=\lambda_{b},\;\bm{q}_{0}\times\bm{n}=\bm{q}_{b}\times\bm{n},\
s_{0}=s_{b},\ \text{ on }\partial T,$

so that $\lambda_{0}\in C(\Omega)$, $s_{0}\in C(\Omega)$ and $\bm{q}_{0}\in
H_{0}(curl;\Omega)$. Hence,

(5.9)
$\displaystyle\nabla\lambda_{0}=\nabla_{w}\lambda_{h},\;\nabla\times\bm{q}_{0}=\nabla_{w}\times\bm{q}_{h}.$

Next, by letting $r=0$ and varying $\bm{v}$ in (5.7) we have

$\varepsilon\nabla_{w}\lambda_{h}+\nabla_{w}\times\bm{q}_{h}=0,$

which, together with (5.9), implies

(5.10) $\varepsilon\nabla\lambda_{0}+\nabla\times\bm{q}_{0}=0.$

From $\bm{q}_{0}\in H_{0}(curl;\Omega)$, we have

$\displaystyle(\varepsilon\nabla\lambda_{0}+\nabla\times\bm{q}_{0},\nabla\lambda_{0})$
$\displaystyle=$
$\displaystyle(\varepsilon\nabla\lambda_{0},\nabla\lambda_{0})+(\nabla\times\bm{q}_{0},\nabla\lambda_{0})$
$\displaystyle=$
$\displaystyle(\varepsilon\nabla\lambda_{0},\nabla\lambda_{0})+\langle\bm{n}\times\bm{q}_{0},\lambda_{0}\rangle$
$\displaystyle=$
$\displaystyle(\varepsilon\nabla\lambda_{0},\nabla\lambda_{0}).$

Thus,

$\displaystyle(\varepsilon\nabla\lambda_{0},\nabla\lambda_{0})=0,$

which gives $\nabla\lambda_{0}=\bm{0}$, and hence $\lambda_{0}\equiv 0$ as a
function with mean value 0. This further leads to $\lambda_{b}\equiv 0$. Thus,
from (5.10) we have

$\nabla\times\bm{q}_{0}=0,\qquad\text{in}\ \Omega.$

Observe that $\bm{q}_{0}$ satisfies

$(\bm{q}_{0},\varepsilon\nabla_{w}r)=0\qquad\forall r\in S_{h},$

which leads to $\bm{q}_{0}\in H(div_{\varepsilon};\Omega)$ and

$\nabla\cdot(\varepsilon\bm{q_{0}})=0,\quad\langle\bm{q}_{0}\cdot{\bm{n}}_{i},1\rangle_{\Gamma_{i}}=0,$

for $i=1,2,\cdots L$. This, together with $\nabla\times\bm{q}_{0}=0$ and
$\bm{q_{0}}\in H_{0}(curl;\Omega)$ shows that $\bm{q}_{0}\equiv 0$, and hence
$\bm{q}_{b}={\bm{n}}\times(\bm{q}_{b}\times{\bm{n}})={\bm{n}}\times 0=0$.

Next, from the Helmholtz decomposition (A.1), we have

$\bm{u}_{h}=\varepsilon^{-1}\nabla\times\bm{\tilde{\psi}}+\nabla\tilde{\phi}+\bm{\tilde{\eta}},$

with some $\bm{\tilde{\eta}}\in\mathbb{H}_{\varepsilon n,0}(\Omega)$ and
$\bm{\tilde{\psi}}\in H_{0}(curl;\Omega)$ satisfying
$\nabla\cdot(\varepsilon\bm{\tilde{\psi}})=0$ and
$\langle\varepsilon\bm{\tilde{\psi}}\cdot{\bm{n}}_{i},1\rangle_{\Gamma_{i}}=0$
for $i=1,\cdots,L$. From $s_{2}(s_{h},s_{h})=0$ we have $s_{0}=s_{b}$ on
${\partial T}$ for each element $T\in{\mathcal{T}}_{h}$ so that $s_{0}\in
H^{1}(\Omega)$. It follows that $\nabla_{w}s_{h}=\nabla s_{0}$. If
$\mathbb{H}_{\varepsilon n,0}(\Omega)$ has dimension 0, then
$\bm{\tilde{\eta}}=0$. In (5.6), by choosing the test function $\phi$ and
$\bm{\psi}$ to be the $L^{2}$ projection of the corresponding function in the
Helmholtz decomposition we arrive at

(5.11)
$\begin{split}0=&(\bm{u}_{h},\varepsilon\nabla_{w}Q_{h}\tilde{\varphi}+\nabla_{w}\times\mathbb{Q}_{h}\bm{\tilde{\psi}})+(\mathbb{Q}_{0}\bm{\tilde{\psi}},\varepsilon\nabla_{w}s_{h})\\\
=&(\bm{u}_{h},{\mathcal{Q}}_{h}\varepsilon\nabla\tilde{\varphi}+{\mathcal{Q}}_{h}\nabla\times\bm{\tilde{\psi}})+(\bm{\tilde{\psi}},\varepsilon\nabla
s_{0})\\\
=&(\bm{u}_{h},\varepsilon\nabla\tilde{\varphi}+\nabla\times\bm{\tilde{\psi}})+(\bm{\tilde{\psi}},\varepsilon\nabla
s_{0})\\\
=&(\varepsilon\bm{u}_{h},\bm{u}_{h}-\bm{\tilde{\eta}})+(\bm{\tilde{\psi}},\varepsilon\nabla
s_{0})\\\
=&(\varepsilon(\bm{u}_{h}-\bm{\tilde{\eta}}),\bm{u}_{h}-\bm{\tilde{\eta}}),\end{split}$

which leads to $\bm{u}_{h}-\bm{\tilde{\eta}}=0$; i.e., $\bm{u}_{h}$ is a
harmonic function. As a harmonic function and piecewise polynomial of degree
$k$, the first term on the left-hand side of (5.6) becomes to be zero for all
test functions $\varphi\in M_{h}$ and $\bm{\psi}\in{\bm{W}}_{h}$. It follows
that $\nabla_{w}s_{h}=0$ so that $\nabla s_{0}=\nabla_{w}s_{h}=0$ and hence
$s_{0}\equiv 0$, so is $s_{b}\equiv 0$. ∎

The following is our main result concerning the solution existence and
uniqueness of the numerical scheme (4.2).

######  Theorem 5.3.

The PDWG finite element scheme (4.2) has one and only one solution for all the
components except $\bm{u}_{h}$. The solution $\bm{u}_{h}$ is unique up to a
discrete harmonic function $\bm{\eta}_{h}\in\mathbb{H}_{\varepsilon
n,0}(\Omega)$ that is a piecewise polynomial of degree $k$.

For the PDWG element of lowest order (i.e., $k=0$), a discrete harmonic
function $\bm{\eta}_{h}\in\mathbb{H}_{\varepsilon n,0}(\Omega)$ would be a
piecewise constant vector field that is continuous across each interior
element interface and has vanishing value on the domain boundary along the
normal direction. It follows that one must have $\bm{\eta}_{h}\equiv 0$, or
equivalently, the PDWG finite element scheme (4.2) has one and only one
solution for all the components.

##  6\. Error Analysis

For the exact solution $\\{{\bm{u}};s=0\\}$ of the div-curl system, we have
from (3.1) and (3.2) that

(6.1)
$\begin{split}B_{h}({\mathcal{Q}}_{h}{\bm{u}},Q_{h}s;\varphi,{\bm{\psi}})=&({\mathcal{Q}}_{h}{\bm{u}},\varepsilon\nabla_{w}\varphi+\nabla_{w}\times{\bm{\psi}})+({\bm{\psi}}_{0},\varepsilon\nabla_{w}Q_{h}s)\\\
=&({\mathcal{Q}}_{h}{\bm{u}},\varepsilon\nabla\varphi_{0}+\nabla\times{\bm{\psi}}_{0})\\\
&+\langle{\mathcal{Q}}_{h}{\bm{u}},\varepsilon{\bm{n}}(\varphi_{b}-\varphi_{0})+({\bm{\psi}}_{0}-{\bm{\psi}}_{b})\times{\bm{n}}\rangle_{\partial{\mathcal{T}}_{h}}\\\
=&({\bm{u}},\varepsilon\nabla\varphi_{0}+\nabla\times{\bm{\psi}}_{0})\\\
&+\langle{\mathcal{Q}}_{h}{\bm{u}},\varepsilon{\bm{n}}(\varphi_{b}-\varphi_{0})+({\bm{\psi}}_{0}-{\bm{\psi}}_{b})\times{\bm{n}}\rangle_{\partial{\mathcal{T}}_{h}}\\\
=&-(\nabla\cdot(\varepsilon{\bm{u}}),\varphi_{0})+(\nabla\times{\bm{u}},{\bm{\psi}}_{0})\\\
&+\langle{\bm{u}},\varepsilon{\bm{n}}(\varphi_{0}-\varphi_{b})+({\bm{\psi}}_{b}-{\bm{\psi}}_{0})\times{\bm{n}}\rangle_{\partial{\mathcal{T}}_{h}}\\\
&+\langle{\mathcal{Q}}_{h}{\bm{u}},\varepsilon{\bm{n}}(\varphi_{b}-\varphi_{0})+({\bm{\psi}}_{0}-{\bm{\psi}}_{b})\times{\bm{n}}\rangle_{\partial{\mathcal{T}}_{h}}\\\
&+\langle{\phi}_{1},\varphi_{b}\rangle_{\partial\Omega}\\\
=&\langle{\phi}_{1},{\varphi}_{b}\rangle_{\partial\Omega}-(f,\varphi_{0})+(\bm{g},\bm{\psi}_{0})\\\
&+\langle{\bm{u}}-{\mathcal{Q}}_{h}{\bm{u}},\varepsilon{\bm{n}}(\varphi_{0}-\varphi_{b})+({\bm{\psi}}_{b}-{\bm{\psi}}_{0})\times{\bm{n}}\rangle_{\partial{\mathcal{T}}_{h}}.\end{split}$

Combining the above equation with the fact that $\lambda=0,\ \bm{q}=0$ we
obtain

(6.2)
$\begin{split}&s_{1}(Q_{h}\lambda-\lambda_{h},\mathbb{Q}_{h}\bm{q}-\bm{q}_{h};\varphi,\bm{\psi})+B({\mathcal{Q}}_{h}{\bm{u}}-{\bm{u}}_{h},Q_{h}s-s_{h};\varphi,{\bm{\psi}})\\\
=&\langle{\bm{u}}-{\mathcal{Q}}_{h}{\bm{u}},\varepsilon{\bm{n}}(\varphi_{0}-\varphi_{b})+({\bm{\psi}}_{b}-{\bm{\psi}}_{0})\times{\bm{n}}\rangle_{\partial{\mathcal{T}}_{h}}.\end{split}$

The second error equation can be easily obtained as follows:

(6.3)
$-s_{2}(Q_{h}s-s_{h},r)+B(\bm{v},r;Q_{h}\lambda-\lambda_{h},\mathbb{Q}_{h}\bm{q}-\bm{q}_{h})=0,$

where we have used the fact that $s=0,\ \bm{q}=0$, and $\lambda=0$.

Denote the error functions by

$e_{{\bm{u}}}={\mathcal{Q}}_{h}{\bm{u}}-{\bm{u}}_{h},\ e_{s}=Q_{h}s-s_{h},\
e_{\lambda}=Q_{h}\lambda-\lambda_{h},\
e_{\bm{q}}=\mathbb{Q}_{h}\bm{q}-\bm{q}_{h}.$

######  Theorem 6.1.

For the numerical solution ${\bm{u}}_{h},\ s_{h},\ \lambda_{h},\ \bm{q}_{h}$
arising from the PDWG scheme (4.2), the following estimate holds true:

(6.4) $\begin{split}{|||}(e_{\lambda},e_{\bm{q}}){|||}+{|||}e_{s}{|||}\lesssim
h^{k+\theta}\|{\bm{u}}\|_{k+\theta},\end{split}$

provided that ${\bm{u}}\in[H^{k+\theta}(\Omega)]^{3}$ for $\theta\in(1/2,1]$.

###### Proof.

From (6.2) and (6.3) we have

$\begin{split}s_{1}(e_{\lambda},e_{\bm{q}};e_{\lambda},e_{\bm{q}})+s_{2}(e_{s},e_{s})=&\langle{\bm{u}}-{\mathcal{Q}}_{h}{\bm{u}},\varepsilon{\bm{n}}(e_{\lambda,0}-e_{\lambda,b})+(e_{\bm{q},b}-e_{\bm{q},0})\times{\bm{n}}\rangle_{\partial{\mathcal{T}}_{h}},\end{split}$

which leads to

(6.5) ${|||}(e_{\lambda},e_{\bm{q}}){|||}^{2}+{|||}e_{s}{|||}^{2}\lesssim
h^{k+\theta}\|{\bm{u}}\|_{k+\theta}{|||}(e_{\lambda},e_{\bm{q}}){|||},$

where $\theta\in(1/2,1]$ and $k$ is the order of polynomials for the finite
element space ${\bm{V}}_{h}$. This completes the proof of the theorem. ∎

To derive an estimate for the error function $e_{\bm{u}}$, we use the
Helmholtz decomposition (A.1) to obtain a $\tilde{\phi}\in H^{1}(\Omega)$,
$\bm{\tilde{\psi}}\in{H}_{0}(curl;\Omega)$, and
$\tilde{\bm{\eta}}\in\mathbb{H}_{\varepsilon n,0}(\Omega)$ such that

(6.6)
$e_{\bm{u}}=\varepsilon^{-1}\nabla\times\bm{\tilde{\psi}}+\nabla\tilde{\phi}+\tilde{\bm{\eta}}.$

Assume the following $H^{\alpha}$-regularity holds true for some fixed
$\alpha\in(1/2,1]$:

(6.7)
$\|\bm{\tilde{\psi}}\|_{\alpha}+\|\tilde{\phi}\|_{\alpha}\lesssim\|e_{\bm{u}}-\tilde{\bm{\eta}}\|_{0}.$

The following is the main convergence result of this paper.

######  Theorem 6.2.

Let ${\bm{u}}$ be a solution of the div-curl system (1.1)-(1.3). Assume that
the Helmholtz decomposition (6.6) has the $H^{\alpha}$-regularity estimate
(6.7). For a numerical solution ${\bm{u}}_{h},\ s_{h},\ \lambda_{h},\
\bm{q}_{h}$ arising from (4.2), there exists a harmonic function
$\bm{\eta}\in\mathbb{H}_{\varepsilon n,0}(\Omega)$ such that the following
estimate holds true:

(6.8)
$\begin{split}\|\varepsilon^{1/2}({\bm{u}}_{h}+\bm{\eta}-{\mathcal{Q}}_{h}{\bm{u}})\|\lesssim
h^{k+\theta+\alpha-1}\|{\bm{u}}\|_{k+\theta},\end{split}$

provided that ${\bm{u}}\in[H^{k+\theta}(\Omega)]^{3}$ for $\theta\in(1/2,1]$.

###### Proof.

From the first error equation (6.2), we have

(6.9)
$\begin{split}s_{1}(e_{\lambda},e_{\bm{q}};\varphi,\bm{\psi})+B(e_{\bm{u}},e_{s};\varphi,{\bm{\psi}})=&\langle{\bm{u}}-{\mathcal{Q}}_{h}{\bm{u}},\varepsilon{\bm{n}}(\varphi_{0}-\varphi_{b})+({\bm{\psi}}_{b}-{\bm{\psi}}_{0})\times{\bm{n}}\rangle_{\partial{\mathcal{T}}_{h}}.\end{split}$

By letting $\varphi=Q_{h}\tilde{\phi}$ and $\bm{\psi}=Q_{h}\bm{\tilde{\psi}}$,
we obtain from Lemma 5.1

(6.10)
$\begin{split}B(e_{\bm{u}},e_{s};\varphi,\bm{\psi})=&(e_{\bm{u}},\varepsilon\nabla_{w}Q_{h}\tilde{\phi}+\nabla_{w}\times\mathbb{Q}_{h}\bm{\tilde{\psi}})+(\mathbb{Q}_{0}\bm{\tilde{\psi}},\varepsilon\nabla_{w}e_{s})\\\
=&(e_{\bm{u}},\varepsilon{\mathcal{Q}}_{h}\nabla_{w}\tilde{\phi}+{\mathcal{Q}}_{h}\nabla_{w}\times\bm{\tilde{\psi}})+(\mathbb{Q}_{0}\bm{\tilde{\psi}},\varepsilon\nabla_{w}e_{s})\\\
=&(e_{\bm{u}},\varepsilon\nabla\tilde{\phi}+\nabla\times\bm{\tilde{\psi}})+(\mathbb{Q}_{0}\bm{\tilde{\psi}},\varepsilon\nabla_{w}e_{s})\\\
=&(\varepsilon
e_{\bm{u}},e_{\bm{u}}-\tilde{\bm{\eta}})+(\varepsilon\mathbb{Q}_{0}\bm{\tilde{\psi}},\nabla_{w}e_{s})\\\
=&(\varepsilon(e_{\bm{u}}-\tilde{\bm{\eta}}),e_{\bm{u}}-\tilde{\bm{\eta}})+(\varepsilon\mathbb{Q}_{0}\bm{\tilde{\psi}},\nabla_{w}e_{s}).\end{split}$

From the definition of the weak gradient we have

$\begin{split}(\varepsilon\mathbb{Q}_{0}\bm{\tilde{\psi}},\nabla_{w}e_{s})=&(\varepsilon\mathbb{Q}_{0}\bm{\tilde{\psi}},\nabla
e_{s,0})+\langle\varepsilon\mathbb{Q}_{0}\bm{\tilde{\psi}}\cdot{\bm{n}},e_{s,b}-e_{s,0}\rangle_{\partial{\mathcal{T}}_{h}}\\\
=&(\varepsilon\bm{\tilde{\psi}},\nabla
e_{s,0})+\langle\varepsilon\mathbb{Q}_{0}\bm{\tilde{\psi}}\cdot{\bm{n}},e_{s,b}-e_{s,0}\rangle_{\partial{\mathcal{T}}_{h}}\\\
=&-(\nabla\cdot(\varepsilon\bm{\tilde{\psi}}),\nabla
e_{s,0})+\langle\varepsilon\bm{\tilde{\psi}}\cdot{\bm{n}},e_{s,0}\rangle_{\partial{\mathcal{T}}_{h}}+\langle\varepsilon\mathbb{Q}_{0}\bm{\tilde{\psi}}\cdot{\bm{n}},e_{s,b}-e_{s,0}\rangle_{\partial{\mathcal{T}}_{h}}\\\
=&\langle\varepsilon\bm{\tilde{\psi}}\cdot{\bm{n}},e_{s,0}-e_{s,b}\rangle_{\partial{\mathcal{T}}_{h}}+\langle\varepsilon\mathbb{Q}_{0}\bm{\tilde{\psi}}\cdot{\bm{n}},e_{s,b}-e_{s,0}\rangle_{\partial{\mathcal{T}}_{h}}\\\
=&\langle\varepsilon(\bm{\tilde{\psi}}-\mathbb{Q}_{0}\bm{\tilde{\psi}})\cdot{\bm{n}},e_{s,0}-e_{s,b}\rangle_{\partial{\mathcal{T}}_{h}}.\end{split}$

Substituting the above into (6.10) then (6.9) yields

$\begin{split}\|\varepsilon^{1/2}(e_{\bm{u}}-\tilde{\bm{\eta}})\|^{2}=&B(e_{\bm{u}},e_{s};\varphi,\bm{\psi})-\langle\varepsilon(\bm{\tilde{\psi}}-\mathbb{Q}_{0}\bm{\tilde{\psi}})\cdot{\bm{n}},e_{s,0}-e_{s,b}\rangle_{\partial{\mathcal{T}}_{h}}\\\
=&\langle{\bm{u}}-{\mathcal{Q}}_{h}{\bm{u}},\varepsilon{\bm{n}}(\varphi_{0}-\varphi_{b})+({\bm{\psi}}_{b}-{\bm{\psi}}_{0})\times{\bm{n}}\rangle_{\partial{\mathcal{T}}_{h}}-s_{1}(e_{\lambda},e_{\bm{q}};\varphi,\bm{\psi})\\\
&-\langle\varepsilon(\bm{\tilde{\psi}}-\mathbb{Q}_{0}\bm{\tilde{\psi}})\cdot{\bm{n}},e_{s,0}-e_{s,b}\rangle_{\partial{\mathcal{T}}_{h}},\end{split}$

which leads to

$\begin{split}\|\varepsilon^{1/2}(e_{\bm{u}}-\tilde{\bm{\eta}})\|^{2}\lesssim&h^{\theta}\|{\bm{u}}-{\mathcal{Q}}_{h}{\bm{u}}\|_{\theta}{|||}(\varphi,{\bm{\psi}}){|||}+{|||}(e_{\lambda},e_{\bm{q}}){|||}{|||}(\varphi,\bm{\psi}){|||}\\\
&+h^{\alpha+(\gamma-1)/2}\|\bm{\tilde{\psi}}\|_{\alpha}{|||}e_{s}{|||}.\end{split}$

It can be proved that

${|||}(\varphi,{\bm{\psi}}){|||}\lesssim
h^{\alpha-1}\|(\tilde{\varphi},\tilde{\bm{\psi}})\|_{\alpha}.$

It follows that

$\begin{split}\|\varepsilon^{1/2}(e_{\bm{u}}-\tilde{\bm{\eta}})\|^{2}\lesssim&h^{\theta}\|{\bm{u}}-Q_{h}{\bm{u}}\|_{\theta}h^{\alpha-1}\|(\tilde{\varphi},\tilde{\bm{\psi}})\|_{\alpha}+{|||}(e_{\lambda},e_{\bm{q}}){|||}h^{\alpha-1}\|(\tilde{\varphi},\tilde{\bm{\psi}})\|_{\alpha}\\\
&+h^{\alpha+(\gamma-1)/2}\|\bm{\tilde{\psi}}\|_{\alpha}{|||}e_{s}{|||}.\end{split}$

so that

$\begin{split}\|\varepsilon^{1/2}(e_{\bm{u}}-\tilde{\bm{\eta}})\|\lesssim&h^{\alpha+\theta-1}\|{\bm{u}}-Q_{h}{\bm{u}}\|_{\theta}+h^{\alpha-1}{|||}(e_{\lambda},e_{\bm{q}}){|||}+h^{\alpha+(\gamma-1)/2}{|||}e_{s}{|||}\\\
\leq&h^{k+\theta+\alpha-1}\|{\bm{u}}\|_{k+\theta},\end{split}$

which gives rise to the error estimate (6.8). ∎

##  7\. Numerical Experiments

In this section, we present some numerical results for the PDWG finite element
method proposed and analyzed in previous sections. For simplicity, we choose
the lowest order PDWG element; i.e., $k=0$ so that the solution $\bm{u}$ is
approximated by discontinuous piecewise constant vector fields. The exact
solution $\bm{u}$ has various regularities ranging from smooth to corner
singular. The computational domain includes convex and non-convex polyhedral
regions; some have cavities or multiple toroidal topology. The implementation
uses an open-source and publicly available MATLAB package iFEM [8]. The
computational domain is first partitioned into cubes, and each cube is further
divided into 6 tetrahedra of equi-volume to form a shape-regular finite
element partition. On each tetrahedral element $T$ with boundary $\partial
T=\cup_{i=1}^{4}F_{i}$, the local finite element space consists of functions
given as follows:

(7.1) $\displaystyle\bm{u}_{h}|_{T}\in\left[P_{0}(T)\right]^{3},$ (7.2)
$\displaystyle
s_{h}|_{T}=\left\\{s_{0},s_{b}\right\\}\in\left\\{P_{0}(T),\Pi_{i=1}^{4}P_{0}(F_{i})\right\\},$
(7.3)
$\displaystyle\varphi_{h}|_{T}=\left\\{\varphi_{0},\varphi_{b}\right\\}\in\left\\{P_{0}(T),\Pi_{i=1}^{4}P_{0}(F_{i})\right\\},$
(7.4)
$\displaystyle\bm{\psi}_{h}|_{T}=\left\\{\bm{\psi}_{0},\bm{\psi}_{b}\right\\}\in\left\\{\left[P_{0}(T)\right]^{3},T_{0}(\partial
T)\right\\}.$

Among the spaces above, $T_{0}(\partial T)$ is tangetial to the boundary and
is given by

$T_{0}(\partial
T):=\\{\bm{\psi}:\bm{\psi}_{F_{i,j}}\in\left[P_{0}(F_{i})\right]^{3}\times\bm{n}_{F_{i}},\;j=1,2\text{
and }i=1,2,3,4\\},$

where $\bm{n}_{F_{i}}$ is the outer unit normal vector to face $F_{i}$. The
basis functions for the first three spaces are straightforward. For
$T_{0}(\partial T)$, on each face $F_{i}$ we choose the normalized vectors
representing the directional vector of any 2 edges ($j=1,2$) among 3 on
$\partial F_{i}$, such that its weak curl is the co-normal vector of this edge
with respect to $F_{i}$ ($i=1,2,3,4$),

$\nabla_{w}\times\left\\{\bm{0},\bm{\psi}_{F_{i,j}}\right\\}=3\bm{\psi}_{F_{i,j}}\times(\nabla\zeta_{i}),$

where $\zeta_{i}$ is the barycentric coordinate associated with the vertex
opposite to face $F_{i}$.

For each test problem, we specify a vector field $\bm{u}$ as the true
solution, while the right-hand side of the div-curl system (1.1)-(1.3) are
computed accordingly. We shall evaluate the following errors for the PDWG
finite element solution:

(7.5)
$\displaystyle\|\varepsilon^{1/2}\mathbf{e}_{\bm{u}}\|:=\|\varepsilon^{1/2}(\bm{u}-\bm{u}_{h})\|,$
(7.6)
$\displaystyle\mathopen{|\mkern-2.5mu|\mkern-2.5mu|}(e_{\lambda},\mathbf{e}_{\bm{q}})\mathclose{|\mkern-2.5mu|\mkern-2.5mu|}:=\big{(}s_{1}(\lambda_{h},\bm{q}_{h};\lambda_{h},\bm{q}_{h})\big{)}^{1/2},$
(7.7)
$\displaystyle\mathopen{|\mkern-2.5mu|\mkern-2.5mu|}e_{s}\mathclose{|\mkern-2.5mu|\mkern-2.5mu|}:=\big{(}s_{2}(s_{h};s_{h})\big{)}^{1/2},$

where the $L^{2}(\Omega)$-norm $\|\varepsilon^{1/2}\mathbf{e}_{\bm{u}}\|$ is
computed by using a higher order Gaussian quadrature on each element.

###  7.1. Example 1

$\Omega=(0,1)^{3}$, $\varepsilon=\operatorname{diag}(3,2,1)$, the true
solution $\bm{u}\in\big{(}H^{1}(\Omega)\big{)}^{3}$ is given by

$\bm{u}(x,y,z)=\left(\begin{array}[]{c}\sin(\pi x)\cos(\pi y)\\\ -\sin(\pi
y)\cos(\pi x)\\\ 0\end{array}\right)+\left(\begin{array}[]{c}x\\\ y\\\
z\end{array}\right).$

The performance of the PDWG finite element solution for this test problem is
illustrated in Table 1. It can be seen that $\mathbf{e}_{\bm{u}}$,
$(e_{\lambda},\mathbf{e}_{\bm{q}})$, and $e_{s}$ all have optimal rate of
convergence for $k=0$ as showen in Theorem 6.2. The plot of the true solution
and the PDWG approximation can be found in Figure 1.

Table 1. Errors and corresponding rates of convergence for Example 1. $1/h$ | $\|\varepsilon^{1/2}\mathbf{e}_{\bm{u}}\|$ | rate | $\mathopen{|\mkern-2.5mu|\mkern-2.5mu|}(e_{\lambda},\mathbf{e}_{\bm{q}})\mathclose{|\mkern-2.5mu|\mkern-2.5mu|}$ | rate | $\mathopen{|\mkern-2.5mu|\mkern-2.5mu|}e_{s}\mathclose{|\mkern-2.5mu|\mkern-2.5mu|}$ | rate
---|---|---|---|---|---|---
2 | 1.64e-1 | – | 2.95e-1 | – | 3.96e-2 | –
4 | 8.16e-2 | 1.01 | 1.68e-1 | 0.82 | 1.86e-2 | 1.09
8 | 3.93e-2 | 1.03 | 8.72e-2 | 0.88 | 8.79e-3 | 1.09
16 | 1.93e-2 | 1.03 | 4.41e-2 | 0.98 | 4.48e-3 | 0.98

(a)

(b)

(c)

Figure 1. The true solution vector field shown in (1(a)) of example 7.1 versus
the WG approximation (1(b)). The vector fields are plotted on four $z=c$
planes. The distribution of $\|\varepsilon^{1/2}\mathbf{e}_{\bm{u}}\|_{T}$
locally is plotted in (1(c)) on the cut plane $x=1/2$ with $h=1/8$.

###  7.2. Example 2

The second example is adopted from [li2018weak] with $\varepsilon=I$ and a
singular solution in $\big{(}H^{1+\frac{2}{3}-\epsilon}(\Omega)\big{)}^{3}$:

$\bm{u}(x,y,z)=\left(\begin{array}[]{c}{x(1-x)}\\\ {y(1-y)}\\\
{r^{2/3}\sin(2\theta)z(1-z)}\end{array}\right),$

in which $r=\sqrt{x^{2}+y^{2}}$, and $\theta=\arctan(y/x)+c$ in the
cylindrical coordinates. Similar to example 1, the result in Table 2 shows
optimal rates of convergence for $\mathbf{e}_{\bm{u}}$, while slightly sub-
optimal in the coarsest two levels for $(e_{\lambda},\mathbf{e}_{\bm{q}})$,
and $e_{s}$, respectively. The plot of the true solution and the PDWG
approximation can be found in Figure 2.

Table 2. Errors and corresponding rates of convergence for Example 2. $1/h$ | $\|\varepsilon^{1/2}\mathbf{e}_{\bm{u}}\|$ | rate | $\mathopen{|\mkern-2.5mu|\mkern-2.5mu|}(e_{\lambda},\mathbf{e}_{\bm{q}})\mathclose{|\mkern-2.5mu|\mkern-2.5mu|}$ | rate | $\mathopen{|\mkern-2.5mu|\mkern-2.5mu|}e_{s}\mathclose{|\mkern-2.5mu|\mkern-2.5mu|}$ | rate
---|---|---|---|---|---|---
2 | 1.13e-1 | – | 2.07e-1 | – | 1.74e-2 | –
4 | 5.20e-2 | 1.12 | 1.20e-1 | 0.78 | 1.05e-2 | 0.73
8 | 2.50e-2 | 1.05 | 6.27e-2 | 0.94 | 5.53e-3 | 0.93
16 | 1.23e-2 | 1.02 | 3.18e-2 | 0.98 | 2.82e-3 | 0.97

(a)

(b)

(c)

Figure 2. The true solution vector field shown in (2(a)) of example 7.2 versus
the PDWG approximation (2(b)). The vector fields are plotted on four $z=c$
planes. The distribution of $\|\varepsilon^{1/2}\mathbf{e}_{\bm{u}}\|_{T}$
locally is plotted in (2(c)) on the cut plane $x=y$ with $h=1/8$.

###  7.3. Example 3

This test problem is defined on
$\Omega=(-1,1)^{2}\times(0,1)\backslash[0,1]\times[-1,0]\times[0,1]$ with
$\varepsilon=I$ and the singular solution in
$\big{(}H^{2/3-\epsilon}(\Omega)\big{)}^{3}$:

$\bm{u}=\nabla\times\left\langle
0,0,r^{2/3}\sin\Big{(}\frac{2}{3}\theta\Big{)}\right\rangle.$

The true solution $\bm{u}$ is a solenoidal vector field (see Figure 3).
Similar to example 2, $r=\sqrt{x^{2}+y^{2}}$ and $\theta=\arctan(y/x)+c$ are
the cylindrical coordinates, where $c$ is chosen such that
$\bm{u}\in{H}(div)\cap H(curl)$. Since $\bm{u}(x,y,z)$ has unbounded
derivatives as $(x,y,z)$ approaches $\\{x=0,y=0\\}\cap\partial\Omega$, the
Gaussian quadrature would yield large error on elements with boundary
intersecting $z$-axis. To overcome this difficulty, we replace the true
solution by its $L^{2}$-projection in the error computation:

$\|\varepsilon^{1/2}\mathbf{e}_{{\mathcal{Q}}_{h}\bm{u}}\|:=\|\varepsilon^{1/2}({\mathcal{Q}}_{h}\bm{u}-\bm{u}_{h})\|.$

It is observed that $\mathbf{e}_{{\mathcal{Q}}_{h}\bm{u}}$,
$(e_{\lambda},\mathbf{e}_{\bm{q}})$, and $e_{s}$ all have maximum possible
rate of convergence ($\approx h^{2/3}$) on the finer meshes (see Table 3). The
plot of the true solution and the PDWG approximation can be found in Figure 4.

(a)

(b)

Figure 3. The true solution vector field shown in (3(a)) of example 7.3 view from above on $z=1/4$ plane together with the level set of its $z$-component. A coarse mesh ($h=1/2$) used in example 7.3 is illustrated in (3(b)).  Table 3. Errors and corresponding rates of convergence for Example 3. $1/h$ | $\|\varepsilon^{1/2}\mathbf{e}_{{\mathcal{Q}}_{h}\bm{u}}\|$ | rate | $\mathopen{|\mkern-2.5mu|\mkern-2.5mu|}(e_{\lambda},\mathbf{e}_{\bm{q}})\mathclose{|\mkern-2.5mu|\mkern-2.5mu|}$ | rate | $\mathopen{|\mkern-2.5mu|\mkern-2.5mu|}e_{s}\mathclose{|\mkern-2.5mu|\mkern-2.5mu|}$ | rate
---|---|---|---|---|---|---
2 | 5.29e-2 | – | 3.35e-1 | – | 4.99e-2 | –
4 | 3.13e-2 | 0.75 | 2.19e-1 | 0.61 | 3.30e-2 | 0.60
8 | 1.91e-2 | 0.72 | 1.41e-1 | 0.63 | 2.14e-2 | 0.62
16 | 1.16e-2 | 0.72 | 9.03e-2 | 0.65 | 1.37e-2 | 0.64

(a)

(b)

(c)

Figure 4. The vector field of ${\mathcal{Q}}_{h}\bm{u}$ is shown in (4(a)) of
example 7.3 versus the PDWG approximation (4(b)). The vector fields are
plotted on several $z=c$ planes. The distribution of
$\|\varepsilon^{1/2}\mathbf{e}_{{\mathcal{Q}}_{h}\bm{u}}\|_{T}$ locally is
plotted in (4(c)) on the cut plane $z=1/4$ when $h=1/8$.

###  7.4. Example 4

This test problem is defined on the domain
$\Omega=(-3/2,1/2)^{3}\backslash[-1,0]^{3}$ such that the domain boundary
$\partial\Omega$ consits of two disjoint connected components
$\Gamma_{0}=\partial(-3/2,1/2)^{3}$ and $\Gamma_{1}=\partial(-1,0)^{3}$. The
true solution is given by

$\bm{u}=\nabla(r^{1/6}),\quad\text{ with }\quad r=\sqrt{x^{2}+y^{2}+z^{2}}.$

It is straightforward to see that $\bm{u}$ is singular near $(0,0,0)$ which is
one of the nonconvex corners of the cavity (see Figure 5), and
$\bm{u}\in[H^{2/3-\epsilon}(\Omega)]^{3}$.

(a)

(b)

Figure 5. (5(a)): plot of the true solution field of example 7.4, view from
above on $z=0$ plane. (5(b)): a coarse mesh ($h=1/2$) used in example 7.4 cut
by the plane $x=0$, with boundary faces of $\Gamma_{1}$ being highlighted in
red.

Table 4 illustrates the convergence of the PDWG method for example 7.4. It can
be seen that $\mathbf{e}_{{\mathcal{Q}}_{h}\bm{u}}$ has the optimal rate of
convergence of $O(h^{2/3})$. On the other hand, the convergence for
$(e_{\lambda},\mathbf{e}_{\bm{q}})$ and $e_{s}$ seems to be approximating the
optimal rate in this numerical test.

For $s=\\{s_{0},s_{b}\\}$, when solving the algebraic system, we seek a
constant $c_{1}$ such that $s_{b}|_{\Gamma_{1}}=c_{1}$ through a simple post-
processing by treating $s_{b}|_{\Gamma_{1}}$ as the fixed DoFs first. Denote
by $\mathbf{U}$ the vector representation of the solution
$(\bm{u}_{h},\lambda_{h},s_{h},\bm{q}_{h})$, and
$\mathbf{U}_{s}=c_{1}\mathbf{S}$, $\mathbf{S}=(0,\cdots,1,\cdots,1,\cdots,0)$,
the indicator vector representing a solution with $s_{b}=1$ on $\Gamma_{1}$
while all other DoFs are zero. Let $A$ be the full stiffness matrix including
from all nodal bases (including boundary faces), while $A^{\rm int}$ be the
stiffness matrix for all the free DoFs: including the interior DoFs for
$\bm{u}_{h}$, $s_{h}$, and $\bm{q}_{h}$, all except 1 fixed DoF for
$\lambda_{h}$. Let $R$ be the restriction operator such that
$R:\mathbf{U}\mapsto\mathbf{U}^{\rm int}$, which is the vector representing
all the aforementioned free DoFs. First we solve the following algebraic
system:

$A^{\rm int}\mathbf{U}^{\rm int}=R\mathbf{F},$

where $\mathbf{F}$ is the full vector of the right hand side. Then the
constant $c_{1}$ is sought by solving the following minimization problem:

$c_{1}=\operatornamewithlimits{argmin}_{a\in\mathbb{R}}\|A(\mathbf{U}^{\rm
int}+a\mathbf{S})-F\|_{\ell^{2}}.$

The plot of the projection of the true solution and the PDWG approximation is
shown in Figure 6.

(a)

(b)

(c)

Figure 6. Example 7.4: The vector field ${\mathcal{Q}}_{h}\bm{u}$ is shown in (6(a)) for Example 7.4 versus the PDWG approximation (6(b)). The vector fields are plotted on several $z=c$ planes. The error distribution of $\|\varepsilon^{1/2}\mathbf{e}_{{\mathcal{Q}}_{h}\bm{u}}\|_{T}$ is plotted in (6(c)) on the cut plane $z=0$ with meshsize $h=1/8$.  Table 4. Errors and corresponding rates of convergence for Example 4. $1/h$ | $\|\varepsilon^{1/2}\mathbf{e}_{{\mathcal{Q}}_{h}\bm{u}}\|$ | rate | $\mathopen{|\mkern-2.5mu|\mkern-2.5mu|}(e_{\lambda},\mathbf{e}_{\bm{q}})\mathclose{|\mkern-2.5mu|\mkern-2.5mu|}$ | rate | $\mathopen{|\mkern-2.5mu|\mkern-2.5mu|}e_{s}\mathclose{|\mkern-2.5mu|\mkern-2.5mu|}$ | rate
---|---|---|---|---|---|---
2 | 1.90-1 | – | 2.51e-1 | – | 2.04e-2 | –
4 | 1.23e-1 | 0.63 | 1.93e-1 | 0.38 | 1.69e-2 | 0.27
8 | 7.78e-2 | 0.66 | 1.35e-1 | 0.51 | 1.24e-2 | 0.44
16 | 4.91e-2 | 0.66 | 9.03e-2 | 0.58 | 8.47e-3 | 0.55

###  7.5. Example 5

In this example we consider singular solutions in the following vector
potential form on a toroidal domain with 1 hole:
$\Omega=\Big{(}(-1,\frac{1}{2})^{2}\backslash[-\frac{1}{2},0]^{2}\times[0,\frac{1}{2}]\Big{)}\times(0,\frac{1}{2})$

$\bm{u}=\nabla\times\left\langle
0,0,r^{\gamma}\sin\big{(}\alpha\theta\big{)}\right\rangle,$

where $r$ and $\theta$ are the cylindrical coordinates defined as in Example
7.3. It can be verified that for $\gamma\neq 1$,
$\bm{u}\in\big{(}H^{\gamma-\epsilon}(\Omega)\big{)}^{3}$, and for $\gamma<1$,
this vector field is singular near a non-convex corner centered at $z$-axis
(see Figure 7).

(a)

(b)

Figure 7. Example 7.5: Plot of the singular true solution vector field with
$\gamma=3/4$ and $\alpha=2$ is shown in (7(a)), view from above on $z=0$
plane. The level set contours are for the $z$-component of the vector
potential. The coarse mesh ($h=1/2$) is shown in (7(b)).

Unlike Example 7.3 in which $\gamma=\alpha$ was set on an L-shaped domain, we
choose $\alpha=2$ in this example, so that the resulting vector field is non-
harmonic for $\gamma\neq\alpha$, and

$\nabla\times\bm{u}=\left\langle
0,0,-\Delta\left(r^{\gamma}\sin\big{(}\alpha\theta\big{)}\right)\right\rangle=\left\langle
0,0,\left(\alpha^{2}-\gamma^{2}\right)\left(x^{2}+y^{2}\right)^{\gamma/2-1}\sin\left(\alpha\theta\right)\right\rangle.$

Consequently, $\nabla\times\bm{u}\notin L^{2}(\Omega)$ if $\gamma\leq 1$.
Nevertheless, due to the unique nature of the PDWG method, we still obtain
noteworthy convergence result for the case of $\gamma\leq 1$. In Table 7, we
have compiled several cases ranging from smooth to singular and plotted the
PDWG approximation vs ${\mathcal{Q}}_{h}\bm{u}$ in Figure 8.

  * •

Regular case $\gamma=1.25$: $\mathbf{e}_{\bm{u}}$ and $e_{s}$ show the optimal
rates of convergence at $O(h)$, while $\mathbf{e}_{{\mathcal{Q}}_{h}\bm{u}}$
shows a superconvergence. $(e_{\lambda},\mathbf{e}_{\bm{q}})$ only shows a
slightly suboptimal rate of convergence during first two refinements, and
optimal thereafter.

  * •

Singular case $\gamma=1$: In this case, we have $\bm{u}\in
H^{1-\epsilon}_{loc}(\Omega)$ and $\bm{u}\notin\bm{H}(\mathbf{curl})$.
$\mathbf{e}_{\bm{u}}$ shows a rate of convergence at $O(h^{0.9})$
asymptotically. Like the smooth case, $\mathbf{e}_{{\mathcal{Q}}_{h}\bm{u}}$
shows a superconvergence with rates higher than $1$.
$(e_{\lambda},\mathbf{e}_{\bm{q}})$ and $e_{s}$ show an optimal rate of
convergence asymptotically.

  * •

Singular case $\gamma=2/3$: In this case, one has $\bm{u}\in
H^{2/3-\epsilon}(\Omega)$, and $\bm{u}\notin\bm{H}(\mathbf{curl})$. Both
$\mathbf{e}_{\bm{u}}$ and $e_{s}$ show optimal rate of convergence at
$O(h^{2/3})$. Similarly to previous two cases,
$\mathbf{e}_{{\mathcal{Q}}_{h}\bm{u}}$ exhibits superconvergence with a rate
of $O(h^{0.9})$. $(e_{\lambda},\mathbf{e}_{\bm{q}})$ shows a convergence with
a slightly suboptimal rate of $O(h^{0.6})$.

Table 5. Errors and the rates of convergence with different $\gamma$ for Example 5 | $1/h$ | $\|\mathbf{e}_{\bm{u}}\|$ | rate | $\|\varepsilon^{1/2}\mathbf{e}_{{\mathcal{Q}}_{h}\bm{u}}\|$ | rate | $\mathopen{|\mkern-2.5mu|\mkern-2.5mu|}(e_{\lambda},\mathbf{e}_{\bm{q}})\mathclose{|\mkern-2.5mu|\mkern-2.5mu|}$ | rate | $\mathopen{|\mkern-2.5mu|\mkern-2.5mu|}e_{s}\mathclose{|\mkern-2.5mu|\mkern-2.5mu|}$ | rate
---|---|---|---|---|---|---|---|---|---
| 2 | 3.96e-1 | – | 1.62e-1 | – | 9.01e-1 | – | 8.23e-2 | –
$\gamma=5/4$ | 4 | 2.09e-1 | 0.92 | 7.69e-2 | 1.07 | 4.94e-1 | 0.87 | 5.12e-2 | 0.68
(smooth) | 8 | 1.07e-1 | 0.96 | 3.23e-2 | 1.25 | 2.63e-1 | 0.91 | 2.70e-2 | 0.93
| 16 | 5.44e-2 | 0.98 | 1.26e-3 | 1.35 | 1.37e-1 | 0.94 | 1.32e-2 | 1.03
| $1/h$ | $\|\mathbf{e}_{\bm{u}}\|$ | rate | $\|\varepsilon^{1/2}\mathbf{e}_{{\mathcal{Q}}_{h}\bm{u}}\|$ | rate | $\mathopen{|\mkern-2.5mu|\mkern-2.5mu|}(e_{\lambda},\mathbf{e}_{\bm{q}})\mathclose{|\mkern-2.5mu|\mkern-2.5mu|}$ | rate | $\mathopen{|\mkern-2.5mu|\mkern-2.5mu|}e_{s}\mathclose{|\mkern-2.5mu|\mkern-2.5mu|}$ | rate
| 2 | 5.34e-1 | – | 2.46e-1 | – | 1.20e0 | – | 1.42e-1 | –
$\gamma=1$ | 4 | 3.06e-1 | 0.80 | 1.28e-1 | 0.94 | 7.11e-1 | 0.76 | 8.95e-2 | 0.66
(singular) | 8 | 1.67e-1 | 0.87 | 5.58e-2 | 1.19 | 4.04e-1 | 0.81 | 4.94e-2 | 0.86
| 16 | 8.82e-2 | 0.88 | 2.55e-2 | 1.13 | 2.25e-1 | 0.85 | 2.56e-2 | 0.95
| $1/h$ | $\|\mathbf{e}_{\bm{u}}\|$ | rate | $\|\varepsilon^{1/2}\mathbf{e}_{{\mathcal{Q}}_{h}\bm{u}}\|$ | rate | $\mathopen{|\mkern-2.5mu|\mkern-2.5mu|}(e_{\lambda},\mathbf{e}_{\bm{q}})\mathclose{|\mkern-2.5mu|\mkern-2.5mu|}$ | rate | $\mathopen{|\mkern-2.5mu|\mkern-2.5mu|}e_{s}\mathclose{|\mkern-2.5mu|\mkern-2.5mu|}$ | rate
| 2 | 8.87e-1 | – | 4.55e-1 | – | 2.05e0 | – | 3.41e-1 | –
$\gamma=2/3$ | 4 | 5.87e-1 | 0.59 | 2.75e-1 | 0.73 | 1.40e0 | 0.55 | 2.39e-1 | 0.51
(singular) | 8 | 3.70e-1 | 0.67 | 1.39e-1 | 0.98 | 9.28e-1 | 0.60 | 1.53e-1 | 0.65
| 16 | 2.34e-1 | 0.66 | 7.51e-2 | 0.89 | 6.01e-1 | 0.62 | 9.51e-2 | 0.68

(a)

(b)

(c)

Figure 8. Example 7.5: The vector field of ${\mathcal{Q}}_{h}\bm{u}$ is shown
in (8(a)) versus the PDWG approximation (8(b)) for $\gamma=2/3$ (singular
case). The vector fields are plotted on several $z=c$ planes. The distribution
of $\|\varepsilon^{1/2}\mathbf{e}_{\bm{u}}\|_{T}$ locally is plotted in (8(c))
on the cut plane $z=1/4$ with meshsize $h=1/8$.

###  7.6. Example 6

In this example we consider a singular solution bearing the same form with
that in Example 7.5 on a toroidal domain with 2 holes:
$\Omega=\Big{[}(-1,\frac{3}{2})^{2}\backslash\big{\\{}[-\frac{1}{2},0]^{2}\times[0,\frac{1}{2}]\cup[\frac{1}{2},1]\times[-\frac{1}{2},0]\big{\\}}\Big{]}\times[0,\frac{1}{2}]$

$\bm{u}=\nabla\times\left\langle
0,0,r_{1}^{\gamma_{1}}\sin\big{(}\alpha\theta_{1}\big{)}+r_{2}^{\gamma_{2}}\sin\big{(}\alpha\theta_{2}\big{)}\right\rangle,$

where $(r_{i},\theta_{i})$ are the cylindrical coordinates centered at a
nonconvex corner of the $i$-th hole, i.e., $r_{1}=\sqrt{x^{2}+y^{2}}$,
$r_{2}=\sqrt{(x-1)^{2}+y^{2}}$, $\theta_{1}=\arctan(y/x)+c_{1}$, and
$\theta_{2}=\arctan\big{(}y/(x-1)\big{)}+c_{2}$. In this example, we choose
$\gamma_{1}=1/2$ and $\gamma_{2}=2/3$ such that the vector field is singular
near the nonconvex corners of both holes (see Figure 9). In fact, the vector
field $\bm{u}$ behaves as $H^{1/2-\epsilon}$-regular in a neighborhood of the
edge $\\{x=0,y=0\\}$, and as $H^{2/3-\epsilon}$-regular in a neighborhood of
the edge $\\{x=1,y=0\\}$. In this example, one has
$\bm{u}\notin\bm{H}(\mathbf{curl})$ as in example 7.5. The convergence results
are shown in Table 6. It can be seen that $\mathbf{e}_{\bm{u}}$,
$(e_{\lambda},\mathbf{e}_{\bm{q}})$, and $e_{s}$ show optimal rates of
convergence at $O(h^{1/2})$. The local error is more prominent near the
nonconvex corner of the first hole where $\bm{u}$ locally is
$H^{1/2-\epsilon}$-regular (see Figure 10). Similarly to previous two cases,
$\mathbf{e}_{{\mathcal{Q}}_{h}\bm{u}}$ shows superconvergence with a rate
approximately at $O(h^{3/4})$.

(a)

(b)

Figure 9. Example 7.6: the singular true solution vector field when $\gamma_{1}=1/2$, $\gamma_{2}=2/3$, and $\alpha=2$ is shown in (9(a)) view from above on $z=0$ plane. The level set contours are for the $z$-component of the vector potential used. The coarse mesh ($h=1/2$) is illustrated in (9(b)). Table 6. Errors and the rates of convergence for Example 6 | $1/h$ | $\|\mathbf{e}_{\bm{u}}\|$ | rate | $\|\varepsilon^{1/2}\mathbf{e}_{{\mathcal{Q}}_{h}\bm{u}}\|$ | rate | $\mathopen{|\mkern-2.5mu|\mkern-2.5mu|}(e_{\lambda},\mathbf{e}_{\bm{q}})\mathclose{|\mkern-2.5mu|\mkern-2.5mu|}$ | rate | $\mathopen{|\mkern-2.5mu|\mkern-2.5mu|}e_{s}\mathclose{|\mkern-2.5mu|\mkern-2.5mu|}$ | rate
---|---|---|---|---|---|---|---|---|---
| 2 | 1.49e0 | – | 8.37e-1 | – | 3.39e0 | – | 6.38e-1 | –
$\gamma_{1}=1/2$ | 4 | 1.04e0 | 0.52 | 5.18e-1 | 0.69 | 2.48e0 | 0.45 | 4.74e-1 | 0.43
$\gamma_{2}=2/3$ | 8 | 6.99e-1 | 0.57 | 2.84e-1 | 0.87 | 1.77e0 | 0.49 | 3.27e-1 | 0.54
| 16 | 4.79e-1 | 0.55 | 1.70e-1 | 0.73 | 1.24e0 | 0.51 | 2.22e-1 | 0.55

(a)

(b)

(c)

Figure 10. Example 7.6: The vector field of ${\mathcal{Q}}_{h}\bm{u}$ is shown
in (10(a)) versus the PDWG approximation (10(b)) for $\gamma_{1}=1/2$ and
$\gamma_{2}=2/3$ (singular case). The vector fields are plotted on several
$z=c$ planes. The distribution of
$\|\varepsilon^{1/2}\mathbf{e}_{\bm{u}}\|_{T}$ locally is plotted in (10(c))
on the cut plane $z=1/4$ with meshsize $h=1/8$.

###  7.7. Example 7

In this example, we report some computational results for a test problem on a
toroidal domain with 1 hole. The numerical data does not show a convergence
due to the presense of a harmonic vector field, according to Theorem 6.2. The
true solution is obtained by combining the ones used in Examples 7.1 and 7.5:

$\bm{u}(x,y,z)=\nabla\times\left\langle
0,0,r^{\gamma}\sin\big{(}\alpha\theta\big{)}\right\rangle+\beta\left(\begin{array}[]{c}\sin(\pi
x)\cos(\pi y)\\\ -\sin(\pi y)\cos(\pi x)\\\ 0\end{array}\right).$

In this test, we choose $\alpha=2$, $\gamma=2/3$, such that
$\bm{u}\in\bm{H}^{2/3-\epsilon}$, while the extra term with coefficient
$\beta$ is smooth thus not affecting the regularity of the solution. Observe
that the extra $\beta$ term is not divergence free. The optimal rate of
convergence should be of order $O(h^{2/3})$ on simply-connected domains
unaffected by the $\beta$ term. As $\beta$ varies, the numerical results do
not demonstrate any convergence for the vector field $\bm{u}$, while optimal
order of convergence are seen for both $(e_{\lambda},\mathbf{e}_{\bm{q}})$ and
$e_{s}$. The numerical performance is in consistency with our theory as
established in Theorem 6.1 for the convergence of
$(e_{\lambda},\mathbf{e}_{\bm{q}})$ and $e_{s}$ and Theorem 6.2 for the
convergence of the vector field $\bm{u}_{h}$ up to a harmonic field. The
vector fields of ${\mathcal{Q}}_{h}\bm{u}$ and $\bm{u}_{h}$ are plotted in
Figure 11 (see (a) and (b)), while their difference
$\bm{\eta}_{h}={\mathcal{Q}}_{h}\bm{u}-\bm{u}_{h}$ is plotted in the same
figure as (c). According to Theorem 6.2, the vector field $\bm{\eta}_{h}$ is
an approximate harmonic field with normal boundary condition.

Table 7. Errors and the rates of convergence with different $\beta$ for Example 7 | $1/h$ | $\|\mathbf{e}_{\bm{u}}\|$ | rate | $\|\varepsilon^{1/2}\mathbf{e}_{{\mathcal{Q}}_{h}\bm{u}}\|$ | rate | $\mathopen{|\mkern-2.5mu|\mkern-2.5mu|}(e_{\lambda},\mathbf{e}_{\bm{q}})\mathclose{|\mkern-2.5mu|\mkern-2.5mu|}$ | rate | $\mathopen{|\mkern-2.5mu|\mkern-2.5mu|}e_{s}\mathclose{|\mkern-2.5mu|\mkern-2.5mu|}$ | rate
---|---|---|---|---|---|---|---|---|---
| 2 | 1.26e0 | – | 8.08e-1 | – | 2.52e0 | – | 3.68e-1 | –
$\gamma=2/3$ | 4 | 8.78e-1 | 0.52 | 6.35e-1 | 0.35 | 1.63e0 | 0.63 | 2.54e-1 | 0.53
$\beta=1$ | 8 | 6.71e-1 | 0.39 | 5.55e-1 | 0.19 | 1.03e0 | 0.66 | 1.58e-1 | 0.68
| 16 | 5.82e-1 | 0.21 | 5.33e-1 | 0.06 | 6.44e-1 | 0.67 | 9.72e-2 | 0.70
| $1/h$ | $\|\mathbf{e}_{\bm{u}}\|$ | rate | $\|\varepsilon^{1/2}\mathbf{e}_{{\mathcal{Q}}_{h}\bm{u}}\|$ | rate | $\mathopen{|\mkern-2.5mu|\mkern-2.5mu|}(e_{\lambda},\mathbf{e}_{\bm{q}})\mathclose{|\mkern-2.5mu|\mkern-2.5mu|}$ | rate | $\mathopen{|\mkern-2.5mu|\mkern-2.5mu|}e_{s}\mathclose{|\mkern-2.5mu|\mkern-2.5mu|}$ | rate
| 2 | 3.42e0 | – | 2.75e0 | – | 5.32e0 | – | 6.52e-1 | –
$\gamma=2/3$ | 4 | 2.89e0 | 0.25 | 2.66e0 | 0.05 | 3.16e0 | 0.75 | 4.05e-1 | 0.69
$\beta=5$ | 8 | 2.71e0 | 0.09 | 2.64e0 | 0.01 | 1.79e0 | 0.82 | 2.21e-1 | 0.87
| 16 | 2.66e0 | 0.03 | 2.64e0 | 0.00 | 1.01e0 | 0.82 | 1.27e-1 | 0.80

(a)

(b)

(c)

Figure 11. Example 7.7: The vector fields of ${\mathcal{Q}}_{h}\bm{u}$ in
(11(a)) and that of the PDWG approximation (11(b)) are visually different when
$\bm{u}$ is not divergence-free on a toroidal domain. The discrete harmonic
field $\bm{\eta}_{h}={\mathcal{Q}}_{h}\bm{u}-\bm{u}_{h}$ is plotted in
(11(c)). The plots are on several $z=c$ planes.

##  Appendix A Helmholtz Decomposition

######  Theorem A.1.

For any vector-valued function ${\bm{u}}\in[L^{2}(\Omega)]^{3}$, there exists
a unique ${\bm{\psi}}\in H_{0}(curl;\Omega),\ \phi\in
H^{1}(\Omega)/\mathbb{R}$, and $\bm{\eta}\in\mathbb{H}_{\varepsilon
n,0}(\Omega)$ such that

(A.1)
$\displaystyle{\bm{u}}=\varepsilon^{-1}\nabla\times{\bm{\psi}}+\nabla\phi+{\bm{\eta}},$
(A.2) $\displaystyle\nabla\cdot(\varepsilon{\bm{\psi}})=0,\
\langle\varepsilon{\bm{\psi}}\cdot{\bm{n}}_{i},1\rangle_{\Gamma_{i}}=0,\
i=1,\dotsc,L.$

Moreover, the following estimate holds true

(A.3) $\|{\bm{\psi}}\|_{H({\rm
curl};\Omega)}+\|\nabla\phi\|_{0}\lesssim(\varepsilon{\bm{u}},{\bm{u}})^{\frac{1}{2}}.$

###### Proof.

The following is a sketch of the proof. Consider the problem of seeking
${\bm{\psi}}\in{\mathbb{W}_{\varepsilon}(\Omega)}$ such that

(A.4)
$(\varepsilon^{-1}\nabla\times{\bm{\psi}},\nabla\times{\bm{\varphi}})=({\bm{u}},\nabla\times{\bm{\varphi}}),\quad\forall\
{\bm{\varphi}}\in{\mathbb{W}_{\varepsilon}(\Omega)}.$

Denote by

$a({\bm{\psi}},{\bm{\varphi}}):=(\varepsilon^{-1}\nabla\times{\bm{\psi}},\nabla\times{\bm{\varphi}})$

the bilinear form defined on ${\mathbb{W}_{\varepsilon}(\Omega)}$. We claim
that $a(\cdot,\cdot)$ is coercive with respect to the $H({\rm
curl};\Omega)$-norm. To this end, it suffices to derive the following estimate

(A.5)
$\|{\bm{v}}\|_{0}\lesssim\|\nabla\times{\bm{v}}\|_{0},\qquad{\bm{v}}\in{\mathbb{W}_{\varepsilon}(\Omega)}.$

In fact, for ${\bm{v}}\in{\mathbb{W}_{\varepsilon}(\Omega)}$, from Theorem 3.4
(Chapter 1) of [12], there exists a vector potential function
${\bm{\omega}}\in[H^{1}(\Omega)]^{3}$ such that

(A.6) $\varepsilon{\bm{v}}=\nabla\times{\bm{\omega}},\
\nabla\cdot{\bm{\omega}}=0,\
\|{\bm{\omega}}\|_{1}\lesssim(\varepsilon{\bm{v}},{\bm{v}})^{\frac{1}{2}}.$

Using the integration by parts and the condition ${\bm{v}}\times{\bm{n}}=0$ on
$\Gamma$, we have

$(\varepsilon{\bm{v}},{\bm{v}})=(\nabla\times{\bm{\omega}},{\bm{v}})=({\bm{\omega}},\nabla\times{\bm{v}}).$

It follows from the Cauchy-Schwarz inequality and (A.6) that

$(\varepsilon{\bm{v}},{\bm{v}})\leq\|{\bm{\omega}}\|_{0}\
\|\nabla\times{\bm{v}}\|_{0}\lesssim(\varepsilon{\bm{v}},{\bm{v}})^{\frac{1}{2}}\
\|\nabla\times{\bm{v}}\|_{0},$

which implies (A.5).

Now from the Lax-Milgram Theorem, there exists a unique
${\bm{\psi}}\in{\mathbb{W}_{\varepsilon}(\Omega)}$ satisfying the equation
(A.4) such that

$\|{\bm{\psi}}\|_{H({\rm
curl};\Omega)}\lesssim(\varepsilon{\bm{u}},{\bm{u}}).$

It is easy to see that ${\mathbb{W}_{\varepsilon}(\Omega)}$ is equivalent to
the following quotient space:

$H_{0}({\rm curl};\Omega)/(\nabla H_{0c}^{1}(\Omega))=\\{{\bm{v}}\in
H_{0}({\rm curl};\Omega):\ (\varepsilon{\bm{v}},\nabla\phi)=0,\ \forall\phi\in
H_{0c}^{1}(\Omega)\\}.$

Thus, by using a Lagrangian multiplier $p\in H_{0c}^{1}(\Omega)$, the problem
(A.4) can be re-formulated as follows: Find ${\bm{\psi}}\in H_{0}({\rm
curl};\Omega)$ and $p\in H_{0c}^{1}(\Omega)$ such that

(A.7)
$\begin{split}(\varepsilon^{-1}\nabla\times{\bm{\psi}},\nabla\times{\bm{\varphi}})+(\varepsilon\nabla
p,{\bm{\varphi}})&=({\bm{u}},\nabla\times{\bm{\varphi}}),\quad\forall\
{\bm{\varphi}}\in H_{0}({\rm curl};\Omega),\\\ ({\bm{\psi}},\varepsilon\nabla
s)&=0,\quad\qquad\qquad\forall\ s\in H_{0c}^{1}(\Omega).\end{split}$

It follows from the first equation of (A.7) that

$\nabla\times({\bm{u}}-\varepsilon^{-1}\nabla\times{\bm{\psi}})-\varepsilon\nabla
p=0.$

Since $p\in H_{0c}^{1}(\Omega)$, the two terms on the left-hand side of the
above equation are orthogonal in the $\varepsilon^{-1}$-weighted
$L^{2}(\Omega)$ norm. Thus, we have

$\nabla\times({\bm{u}}-\varepsilon^{-1}\nabla\times{\bm{\psi}})=0,$

which gives

${\bm{u}}-\varepsilon^{-1}\nabla\times{\bm{\psi}}\in H^{0}({\rm
curl};\Omega).$

Thus, there exist unique $\phi\in H^{1}(\Omega)/\mathbb{R}$ and
$\bm{\eta}\in\mathbb{H}_{\varepsilon n,0}(\Omega)$ such that

${\bm{u}}-\varepsilon^{-1}\nabla\times{\bm{\psi}}=\nabla\phi+\bm{\eta},$

which completes the proof of the theorem. ∎

## Acknowledgement

We would like to express our gratitude to Dr. Long Chen (UCI) for his valuable
discussion and suggestions.

## References

  * [1] G. Auchmuty and J.C. Alexander, $L^{2}$ Well-posedness of planar div-curl systems, Archive for Rational Mechanics and Analysis, 20 (2001), 160(2), pp. 91-134.
  * [2] I. Babus̆ka, The finite element method with penalty, Math. Comp., 27 (1973), 221-228.
  * [3] R. Bensow and M.G. Larson, Discontinuous least-squares finite element method for the div-curl problem, Numer. Math. 101 (2005), pp. 601-617.
  * [4] P.B. Bochev, K. Peterson and C.M. Siefert, Analysis and computation of compatible least-squares methods for div-curl equations, SIAM J. Numer. Anal., 49 (2011), pp. 159-181.
  * [5] J. Bramble and J.E. Pasciak, A new approximation technique for div-curl systems, Math. Comp. 73 (2004), pp. 1739-1762.
  * [6] F. Brezzi, On the existence, uniqueness, and approximation of saddle point problems arising from Lagrange multipliers, RAIRO, 8 (1974), 129-151.
  * [7] W. Cao and C. Wang, New primal-dual weak Galerkin finite element methods for convection-diffusion problems, arxiv: submit/3011322
  * [8] L. Chen, iFEM: an innovative finite element methods package in MATLAB, UC Irvine, 2009, https://github.com/lyc102/ifem.
  * [9] P.G. Ciarlet, The Finite Element Method for Elliptic Problems, North-Holland, 1978.
  * [10] D.M. Copeland, J. Gopalakrishnan and J.E. Pasciak, A mixed method for axisymmetric div-curl systems, Math. Comput., 77 (2008), pp. 1941-1965.
  * [11] K.O. Friedrichs, Differential forms on Riemannian manifolds, Comm. Pure & Applied Math., vol. 8(1955), 551-590. MR0087763 (19:407a).
  * [12] V. Girault and P-A Raviart, Finite Element Methods for Navier-Stokes Equations: Theory and Algorithms, Springer-Verlag Berlin Heidelberg, 1986.
  * [13] D. Li, C. Wang and J. Wang, Primal-dual weak Galerkin finite element methods for linear convection equations in non-divergence form, arXiv: 1910.14073.
  * [14] J. Li and Y. Huang, Time-domain finite element methods for Maxwell’s equations in metamaterials, Springer, 2013.
  * [15] J. Li, X. Ye and S. Zhang, A weak Galerkin least-squares finite element method for div-curl systems, Journal of Computational Physics, Vol. 363, pp. 79-86. 2018.
  * [16] K. Lipnikov, G. Manzini, F. Brezzi and A. Buffa, The mimetic finite difference method for the 3D magnetostatic field problems on polyhedral meshes, J. Comput. Phys., 230 (2011), pp. 305-328.
  * [17] Y. Liu and J. Wang, A primal-dual weak Galerkin method for div-curl systems with low-regularity solutions, arXiv:2003.11795v2.
  * [18] R. Nicolaides and X. Wu, Covolume solutions of three-dimensional div-curl equations, SIAM J. Numer. Anal. 34 (1997) 2195-2203.
  * [19] J. Saranen, On generalized harmonic fields in domains with anisotropic nonhomogeneous media, J. Math. Anal. Appl., 88 (1982), pp. 104-115.
  * [20] C. Wang, A new primal-dual weak Galerkin finite element method for ill-posed elliptic Cauchy problems, Journal of Computational and Applied Mathematics, 2019, available online.
  * [21] C. Wang and J. Wang, Discretization of div-curl systems by weak Galerkin finite element methods on polyhedral partitions, J. Sci. Comput. 68 (2016) 1144-1171.
  * [22] C. Wang and J. Wang, A primal-dual finite element method for first-order transport problems, Journal of Computational Physics, Vol. 417, 109571, 2020.
  * [23] C. Wang and J. Wang, A primal-dual weak Galerkin finite element method for second order elliptic equations in non-divergence form, Math. Comp., 87 (2018) 515-545.
  * [24] C. Wang and J. Wang, A primal-dual weak Galerkin finite element method for Fokker-Planck type equations, SIAM J. Numer. Anal., 58(5), 2632-2661. 2020. arXiv:1704.05606.
  * [25] C. Wang and J. Wang, Primal-dual weak Galerkin finite element methods for elliptic Cauchy problems, Computers and Mathematics with Applications, vol 79(3), pp. 746-763, 2020.
  * [26] C. Wang and L. Zikatanov, Low regularity primal-dual weak Galerkin finite element methods for convection-diffusion equations, arXiv:1901.06743.
  * [27] J. Wang and X. Ye, A weak Galerkin mixed finite element method for second-order elliptic problems, Math. Comp., vol. 83, pp. 2101-2126, 2014.