New Conforming Finite element divdiv Complexes in Three Dimensions

This paper considers the construction of conforming finite element sub-complexes of the following divdiv complex on $\Omega\subseteq\mathbb{R}^{3}$, see [15, 2]

(1.1)		$\displaystyle RT\stackrel{{\scriptstyle\subseteq}}{{\rightarrow}}H^{1}(\Omega;\mathbb{R}^{3})\stackrel{{\scriptstyle\operatorname{dev}\operatorname{grad}}}{{\longrightarrow}}H(\operatorname{sym}\operatorname{curl},\Omega;$	$\displaystyle\mathbb{T})\stackrel{{\scriptstyle\operatorname{sym}\operatorname{curl}}}{{\longrightarrow}}H(\operatorname{div}\operatorname{div},\Omega;\mathbb{S})$
			$\displaystyle\stackrel{{\scriptstyle\operatorname{div}\operatorname{div}}}{{\longrightarrow}}L^{2}(\Omega;\mathbb{R})\stackrel{{\scriptstyle}}{{\rightarrow}}0,$

where $RT:=\{a\mathbf{x}+\bm{b}:\,a\in\mathbb{R},\bm{b}\in\mathbb{R}^{3}\}$, the spaces $H^{1}(\Omega;\mathbb{R}^{3})$ and $L^{2}(\Omega;\mathbb{R})$ are the vector-valued and scalar-valued Sobolev spaces, respectively, and

	$\displaystyle H(\operatorname{sym}\operatorname{curl},\Omega;\mathbb{T}):=$	$\displaystyle\left\{\bm{\sigma}\in L^{2}(\Omega;\mathbb{T}):\operatorname{sym}\operatorname{curl}\bm{\sigma}\in L^{2}(\Omega;\mathbb{S})\right\},$
	$\displaystyle H(\operatorname{div}\operatorname{div},\Omega;\mathbb{S}):=$	$\displaystyle\left\{\bm{\sigma}\in L^{2}(\Omega;\mathbb{S}):\operatorname{div}\operatorname{div}\bm{\sigma}\in L^{2}(\Omega;\mathbb{R})\right\}.$

Here $\mathbb{T}$ and $\mathbb{S}$ denote the spaces of traceless and symmetric matrices in three dimensions, respectively, the operator $\operatorname{sym}\operatorname{curl}$ is the symmetric part of the row-wise $\operatorname{curl}$ operator, and the inner divergence operator of $\operatorname{div}\operatorname{div}$ is applied by row resulting a column vector for which the outer divergence operator is applied. Such a complex is exact provided that the domain $\Omega$ is contractible and Lipschitz [15, 2]. Particularly, conforming finite element spaces $\Sigma_{h}\subseteq H(\operatorname{div}\operatorname{div},\Omega;\mathbb{S})$, $V_{h}\subseteq H^{1}(\Omega;\mathbb{R}^{3})$, $\mathcal{U}_{h}\subseteq H(\operatorname{sym}\operatorname{curl},\Omega;\mathbb{T})$, and $Q_{h}\subseteq L^{2}(\Omega;\mathbb{R})$ are constructed, such that the following discrete divdiv complex

(1.2)

$$RT\stackrel{{\scriptstyle\subseteq}}{{\longrightarrow}}V_{h}\stackrel{{\scriptstyle\operatorname{dev}\operatorname{grad}}}{{\longrightarrow}}\mathcal{U}_{h}\stackrel{{\scriptstyle\text{ sym curl }}}{{\longrightarrow}}\Sigma_{h}\stackrel{{\scriptstyle\text{ div div }}}{{\longrightarrow}}Q_{h}\rightarrow 0$$

is exact.

One of the applications of the finite element divdiv complex of (1.2) is to solve the fourth order problems [5, 7, 12, 15]. Besides, the construction of the finite element divdiv complex of (1.2) is closely related to the mixed formulation of the so-called linearized Einstein-Bianchi system [16]. The first finite element sub-complex of the gradgrad complex (the dual complex of the divdiv complex) was constructed in [10], and the finite element spaces were employed to solve the linearized Einstein-Bianchi system within the mixed form. Recently, two finite element divdiv complexes of (1.2) on tetrahedral grids were constructed in [11, 7], and the associated finite element spaces can be used to discretize the linearized Einstein-Bianchi system within the dual formulation introduced in [11]. Both used the $H(\operatorname{div}\operatorname{div},\Omega;\mathbb{S})$ conforming finite element on tetrahedral grids from [7]. However, it is arduous to compute the basis functions of the finite element spaces in [11, 7], since the degrees of freedom(DOFs) of those finite elements are in some sense complicated. More precisely, the design of the $H(\operatorname{div}\operatorname{div},\Omega;\mathbb{S})$ conforming finite element in [7] was based on Green’s identity [7, Lemma 4.1]

(1.3)

$\displaystyle\begin{split}(\operatorname{div}\operatorname{div}\bm{\sigma},q)_{K}=&(\bm{\sigma},\nabla^{2}q)_{K}-\sum_{f\in\mathscr{F}(K)}\sum_{e\in\mathscr{E}(f)}(\bm{n}_{f,e}^{\intercal}\bm{\sigma}\bm{n}_{f},q)_{f}\\ &+\sum_{f\in\mathscr{F}(K)}[(\bm{n}_{f}^{\intercal}\bm{\sigma}\bm{n}_{f},\frac{\partial q}{\partial\bm{n}_{f}})_{f}-(2\operatorname{div}_{f}(\bm{\sigma}\bm{n}_{f})+\frac{\partial(\bm{n}_{f}^{\intercal}\bm{\sigma}\bm{n}_{f})}{\partial\bm{n}_{f}},q)_{f}]\end{split}$

on tetrahedron $K$. Here $\bm{\sigma}$ is a symmetric matrix valued function, $q$ is a scalar function, $\mathscr{F}(K)$ is the set of all faces of the tetrahedron $K$, $\mathscr{E}(f)$ is the set of all edges of face $f$, $\bm{n}_{f}$ is the unit normal vector of $f$, $\bm{n}_{f,e}$ denotes the unit normal vector of $e$ that parallels to $f$, and $\operatorname{div}_{f}$ is the surficial divergence. Based on (1.3), besides the continuity of $\bm{n}_{f,e}^{\intercal}\bm{\sigma}\bm{n}_{f}$ across each edge $e$ and $\bm{n}_{f}^{\intercal}\bm{\sigma}\bm{n}_{f}$ across each interface $f$, the continuity of $2\operatorname{div}_{f}(\bm{\sigma}\bm{n}_{f})+\frac{\partial(\bm{n}_{f}^{\intercal}\bm{\sigma}\bm{n}_{f})}{\partial\bm{n}_{f}}$ across each interface $f$ is imposed on the functions of the $H(\operatorname{div}\operatorname{div},\Omega;\mathbb{S})$ conforming finite element space in [7]. As a result, the last continuity brings a set of DOFs which combine two terms $2\operatorname{div}_{f}(\bm{\sigma}\bm{n}_{f})$ and $\frac{\partial(\bm{n}_{f}^{\intercal}\bm{\sigma}\bm{n}_{f})}{\partial\bm{n}_{f}}$ on the boundary of $K$. Such a combination makes it difficult to figure out the explicit basis functions of the corresponding $H(\operatorname{div}\operatorname{div},\Omega;\mathbb{S})$ conforming finite element space. In addition, such a combination is inherited by both the $H^{1}(\Omega;\mathbb{R}^{3})$ and $H(\operatorname{sym}\operatorname{curl},\Omega;\mathbb{T})$ conforming finite element spaces in the finite element complex. In a recent work [12], new $H(\operatorname{div}\operatorname{div},\Omega;\mathbb{S})$ conforming finite element spaces with explicit basis functions were constructed on both triangular and tetrahedral grids in a unified way. Instead of (1.3), the design of [12] was based on the following identity

$$\begin{split}(\operatorname{div}\operatorname{div}\bm{\sigma},q)_{K}=(\bm{\sigma},\nabla^{2}q)_{K}-\sum_{f\in\mathscr{F}(K)}(\bm{\sigma}\bm{n}_{f},\nabla q)_{f}+\sum_{f\in\mathscr{F}(K)}(\bm{n}_{f}^{\intercal}\operatorname{div}\bm{\sigma},q)_{f}.\end{split}$$

The continuity of $\bm{\sigma}\bm{n}_{f}$ and $\bm{n}_{f}^{\intercal}\operatorname{div}\bm{\sigma}$ across each interface $f$ was imposed on the functions of the $H(\operatorname{div}\operatorname{div},\Omega;\mathbb{S})$ conforming finite element space in [12]. Consequently, the element follows by strengthening the continuity of the normal component of the divergences of the functions of the $H(\operatorname{div},\Omega;\mathbb{S})$ conforming finite elements proposed in [14, 9, 13]. Here $H(\operatorname{div},\Omega;\mathbb{S})$ is a space of square-integrable tensors with square-integrable divergence, taking values in the space $\mathbb{S}$ of symmetric matrices. As it can be seen from [12] that because of a crucial structure of the $H(\operatorname{div}\operatorname{div},\Omega;\mathbb{S})$ conforming finite element space, inherited from that of the $H(\operatorname{div},\Omega;\mathbb{S})$ finite element space, the basis functions can be easily written out. In addition, the $H(\operatorname{div}\operatorname{div},\Omega;\mathbb{S})$ conforming finite element space in [12] is a little bit smoother than that in [7]. Recently, the construction in [12] has been generalized to arbitrary dimension [6].

This paper is devoted to construct conforming finite element complexes of (1.2) on both cuboid and tetrahedral grids in three dimensions. The key part is to characterize the continuity of the finite element spaces appearing in the discrete complexes of (1.2). Motivated by the sufficient continuity condition proposed in [12], which is actually a characterization of the continuity of $H(\operatorname{div}\operatorname{div},\Omega;\mathbb{S})\cap H(\operatorname{div},\Omega;\mathbb{S})$, a new $H(\operatorname{div}\operatorname{div},\Omega;\mathbb{S})$ conforming finite element space $\Sigma_{h}$ is constructed in this paper. The enhanced smoothness of $\Sigma_{h}$ brings corresponding enhancement of regularity for the finite element subspaces $\mathcal{U}_{h}\subseteq H(\operatorname{sym}\operatorname{curl},\Omega;\mathbb{T})$ and $V_{h}\subseteq H^{1}(\Omega;\mathbb{R}^{3})$. In fact, a global regularity $\mathcal{U}_{h}\subseteq H^{*}(\operatorname{sym}\operatorname{curl},\Omega;\mathbb{T})$ is imposed for the $H(\operatorname{sym}\operatorname{curl},\Omega;\mathbb{T})$ conforming finite element space, where $H^{*}(\operatorname{sym}\operatorname{curl},\Omega;\mathbb{T})$ is a subspace of $H(\operatorname{sym}\operatorname{curl},\Omega;\mathbb{T})$ such that the column-wise divergences of matrix valued functions belong to the space $H(\operatorname{curl},\Omega;\mathbb{R}^{3})$. Consequently, a global regularity $V_{h}\subseteq H^{1}(\operatorname{div},\Omega;\mathbb{R}^{3})$ is required for the $H^{1}(\Omega;\mathbb{R}^{3})$ conforming finite element space, where $H^{1}(\operatorname{div},\Omega;\mathbb{R}^{3})$ is a subspace of $H^{1}(\Omega;\mathbb{R}^{3})$ such that the divergences of vector valued functions are equally in the space $H^{1}(\Omega;\mathbb{R})$. As a consequence, the conforming finite element divdiv complexes of (1.2) on cuboid grids and tetrahedral grids in three dimensions constructed in this paper are sub-complexes of the divdiv complex with the following enhanced smoothness:

(1.4)		$\displaystyle RT\stackrel{{\scriptstyle\subseteq}}{{\rightarrow}}H^{1}(\operatorname{div},\Omega;$	$\displaystyle\mathbb{R}^{3})\stackrel{{\scriptstyle\operatorname{dev}\operatorname{grad}}}{{\longrightarrow}}H^{\ast}(\operatorname{sym}\operatorname{curl},\Omega;\mathbb{T})$
			$\displaystyle\stackrel{{\scriptstyle\operatorname{sym}\operatorname{curl}}}{{\longrightarrow}}H(\operatorname{div}\operatorname{div},\Omega;\mathbb{S})\cap H(\operatorname{div},\Omega;\mathbb{S})\stackrel{{\scriptstyle\operatorname{div}\operatorname{div}}}{{\longrightarrow}}L^{2}(\Omega;\mathbb{R})\stackrel{{\scriptstyle}}{{\rightarrow}}0.$

On cuboid grids, the spaces $Q_{h}$, $\Sigma_{h}$, $\mathcal{U}_{h}$, and $V_{h}$ are noted as $Q_{k-2,\Box}$, $\Sigma_{k,\Box}$, $\mathcal{U}_{k,\Box}$, and $V_{k,\Box}$ with $k\geq 3$, respectively, where $Q_{k-2,\Box}$ is the space of discontinuous piecewise polynomials of degree no more than $k-2$ for each variable. For the spaces $\Sigma_{k,\Box}$, $\mathcal{U}_{k,\Box}$, and $V_{k,\Box}$, the associated shape function spaces will be determined by the polynomial de Rham complex and the polynomial divdiv complex, and the corresponding DOFs are identified by imposing the associated continuity requirements sufficiently. On tetrahedral grids, the spaces $Q_{h}$, $\Sigma_{h}$, $\mathcal{U}_{h}$, and $V_{h}$ are noted as $Q_{k-2,\triangle}$, $\Sigma_{k,\triangle}$, $\mathcal{U}_{k+1,\triangle}$, and $V_{k+2,\triangle}$ with $k\geq 4$, respectively, where $Q_{k-2,\triangle}$ is the space of discontinuous piecewise polynomials of degree no more than $k-2$, the space $\Sigma_{k,\triangle}$ is a modification of the $H(\operatorname{div}\operatorname{div},\Omega;\mathbb{S})$ conforming finite element in [12], the spaces $\mathcal{U}_{k+1,\triangle}$ and $V_{k+2,\triangle}$ are newly proposed herein.

The rest of the paper is organized as follows. Section 2 introduces the notation. Section 3 constructs the finite element divdiv complex on cuboid grids. Section 4 constructs the finite element divdiv complex on tetrahedral grids. The exactness of the discrete complexes will be proved in Section 5.

Let $\Omega$ be a contractible domain with Lipschitz boundary $\partial\Omega$ of $\mathbb{R}^{3}$.  Let $\{\mathcal{T}_{h}\}$ denote a family of shape regular cuboid or tetrahedral grids on $\Omega$ with the mesh size $h$.  Let $\mathscr{V}$, $\mathscr{E}$, and $\mathscr{F}$ denote the set of vertices, edges, and faces, respectively.  Given element $K\in\mathcal{T}_{h}$, let $\mathscr{V}(K)$, $\mathscr{E}(K)$, and $\mathscr{F}(K)$ denote the sets of vertices, edges, and faces of $K$, respectively. Given an edge $e\in\mathscr{E}(K)$, the unit tangential vector $\bm{t}_{e}$, as well as two unit normal vectors $\bm{n}^{\pm}_{e}$ are fixed. For a facet $f\in\mathscr{F}(K)$, the unit normal vector ${\bm{n}_{f}}$ as well as two unit tangential vectors $\bm{t}^{\pm}_{f}$ are fixed. The jump of $u$ across an interface $f\in{\mathscr{F}}$ shared by neighboring elements $K^{\pm}$ is defined by

$$\left[u\right]_{f}:=u|_{K^{+}}-u|_{K^{-}}.$$

When it comes to any boundary face $f\subseteq\partial\Omega$, the jump $[\cdot]_{f}$ reduces to the trace. Let $\nabla_{pw}^{k}$ denote the elementwise derivative of order $k$ for piecewise functions.

Denote by $\mathbb{M}$ the space of $3\times 3$ real matrices, and let $\mathbb{S}$, $\mathbb{K}$, and $\mathbb{T}$ be the subspaces of symmetric, skew-symmetric, and traceless matrices, respectively. Let $\bm{I}$ denote the $3\times 3$ identity, and $\bm{I}_{1}:=(1,0,0)^{\intercal}$, $\bm{I}_{2}:=(0,1,0)^{\intercal}$, as well as $\bm{I}_{3}:=(0,0,1)^{\intercal}$. For matrix valued function $\bm{\sigma}:\Omega\rightarrow\mathbb{M}$, define the symmetric part and traceless part as follows, respectively,

$$\operatorname{sym}\bm{\sigma}:=\frac{1}{2}(\bm{\sigma}+\bm{\sigma}^{\intercal}),\quad\operatorname{dev}\bm{\sigma}:=\bm{\sigma}-\frac{1}{3}\operatorname{tr}(\bm{\sigma})\bm{I}.$$

For simplicity of presentation, for $\mathbf{x}=(x,y,z)^{\intercal}\in\Omega$, let $x_{i}$ be $x,y,z$ when $i$ takes $1,2$ and $3$, respectively, and $\partial/\partial{x_{i}}$ be abbreviated as $\partial_{i}$. For a vector field $\bm{v}=(v_{1},v_{2},v_{3})^{\intercal}$, the gradient applies by row to produce a matrix valued function, namely $\operatorname{grad}\bm{v}:=\nabla\bm{v}:=(\partial_{j}v_{i})_{3\times 3}$. Define the symmetric gradient by $\varepsilon(\bm{v}):=\operatorname{sym}(\nabla\bm{v})$ and define a skew symmetric matrix

$$\operatorname{mspn}\bm{v}:=\left(\begin{array}[]{ccc}0&-v_{3}&v_{2}\\ v_{3}&0&-v_{1}\\ -v_{2}&v_{1}&0\end{array}\right).$$

Given a facet $f\in\mathscr{F}$ with the unit normal vector $\bm{n}_{f}$, for a vector $\bm{v}\in\mathbb{R}^{3}$ as well as a scalar function $q$, define

$$\Pi_{f}\bm{v}:=(\bm{n}_{f}\times\bm{v})\times\bm{n}_{f},\quad\nabla_{f}q:=\Pi_{f}\nabla q.$$

Define the surface symmetric gradient $\varepsilon_{f}$ by

$$\varepsilon_{f}(\bm{v}):=\operatorname{sym}(\nabla_{f}(\Pi_{f}\bm{v})),$$

and the surface rot operator is defined by

$$\operatorname{rot}_{f}\bm{v}:=(\bm{n}_{f}\times\nabla)\cdot\bm{v}=\bm{n}_{f}\cdot\operatorname{curl}\bm{v}.$$

For a matrix field, the operators $\operatorname{curl}$ and $\operatorname{div}$ apply by row, and $\operatorname{curl}^{*}$ and $\operatorname{div}^{*}$ apply by column to produce a matrix field and vector field, respectively. Given a plane $f$ with unit normal vector $\bm{n}_{f}$, denote

(2.1)

$\displaystyle\mathcal{Q}_{f}:=\bm{I}-\bm{n}_{f}\bm{n}_{f}^{\intercal}$

the orthogonal projection onto $f$.

Given a matrix valued function $\bm{\sigma}$, define the symmetric projection $\mathcal{Q}_{f,\operatorname{sym}}\bm{\sigma}:=\operatorname{sym}(\mathcal{Q}_{f}\bm{\sigma})$, and for symmetric $\bm{\sigma}$, define $\Lambda_{f}(\bm{\sigma}):f\rightarrow\mathcal{Q}_{f}\mathbb{S}\mathcal{Q}_{f}$ by

(2.2)

$\displaystyle\Lambda_{f}(\bm{\sigma}):=\mathcal{Q}_{f}(2\varepsilon(\bm{\sigma}\bm{n}_{f})-\frac{\partial\bm{\sigma}}{\partial\bm{n}_{f}})\mathcal{Q}_{f}.$

Assume that $D\subseteq\mathbb{R}^{3}$ is a bounded and topologically trivial domain. Let standard notation $H^{m}(D;X)$ denote the Sobolev space consisting of functions within domain $D$, taking values in space $X$, and with all derivatives of order at most $m$ square-integrable. In the case $m=0$, set $H^{0}(D;X)=L^{2}(D;X)$. The $L^{2}$ scalar product over $D$ is denoted as $(\cdot,\cdot)_{D}$, $\|{\cdot}\|_{0,D}$ denotes the $L^{2}$ norm over a set $D$, and $\|{\cdot}\|_{0}$ abbreviates $\|{\cdot}\|_{0,\Omega}$. Similarly, let $C^{m}(D;X)$ denote the space of $m$-times continuously differentiable functions, taking values in $X$. In this paper, $X$ could be $\mathbb{M},\mathbb{S},\mathbb{T},\mathbb{R}$ or $\mathbb{R}^{3}$. Define

	$\displaystyle H(\operatorname{sym}\operatorname{curl},D;\mathbb{T}):=$	$\displaystyle\left\{\bm{\sigma}\in L^{2}(D;\mathbb{T}):\operatorname{sym}\operatorname{curl}\bm{\sigma}\in L^{2}(D;\mathbb{S})\right\},$
	$\displaystyle H(\operatorname{div}\operatorname{div},D;\mathbb{S}):=$	$\displaystyle\left\{\bm{\sigma}\in L^{2}(D;\mathbb{S}):\operatorname{div}\operatorname{div}\bm{\sigma}\in L^{2}(D;\mathbb{R})\right\}.$

Let $P_{k}(D;X)$ denote the space of polynomials of degree no more than $k$ on $D$, taking values in the space $X$. Let $Q_{k}(D;X)$ denote the space of polynomials on $D$ of degree no more than $k$ for each variable, taking values in the space $X$. It will be convenient to define

$$P_{k_{1}}(x_{i})\cdot P_{k_{2}}(x_{j})\cdot P_{k_{3}}(x_{l}):=\left\{\sum_{0\leq s_{1}\leq k_{1},0\leq s_{2}\leq k_{2}\atop 0\leq s_{3}\leq k_{3}}c_{s_{1}s_{2}s_{3}}\,x_{i}^{s_{1}}x_{j}^{s_{2}}x_{l}^{s_{3}}:\,c_{s_{1}s_{2}s_{3}}\in\mathbb{R}\right\}.$$

Here $s_{1}$, $s_{2}$, $s_{3}$, $k_{1}$, $k_{2}$, and $k_{3}$ are nonnegative integers, and $\{i,j,l\}$ is a permutation of $\{1,2,3\}$. In particular, for $\mathbf{x}=(x_{1},x_{2},x_{3})^{\intercal}\in D$, denote

$$P_{k_{1}k_{2}k_{3}}(D):=P_{k_{1}}(x_{1})\cdot P_{k_{2}}(x_{2})\cdot P_{k_{3}}(x_{3}).$$

Given two spaces $\mathcal{Q}_{1}(D)$ and $\mathcal{Q}_{2}(D)$ with $\mathcal{Q}_{2}(D)\subseteq\mathcal{Q}_{1}(D)$, let $\mathcal{Q}_{1}(D)/\mathcal{Q}_{2}(D)$ denote the space of $q\in\mathcal{Q}_{1}(D)$ such that $(q,p)_{D}=0$ for any $p\in\mathcal{Q}_{2}(D)$.

For $\mathbf{x}=(x,y,z)^{\intercal}\in D$, denote

(2.3)

$\displaystyle RM(D):=\left\{\begin{pmatrix}c_{1}-c_{4}y-c_{5}z\\ c_{2}+c_{4}x-c_{6}z\\ c_{3}+c_{5}x+c_{6}y\\ \end{pmatrix}:c_{1},c_{2},c_{3},c_{4},c_{5},c_{6}\in\mathbb{R}\right\}$

as the 6-dimensional space of the rigid motions. Besides,

$$RT(D):=\{a\mathbf{x}+\bm{b}:\,a\in\mathbb{R},\bm{b}\in\mathbb{R}^{3}\}$$

with ${\rm{dim}}\,RT(D)=4$ is the shape function space of the lowest order Raviart-Thomas element[4]. If the domain $D$ is clear, the spaces $RM(D)$ and $RT(D)$ are recorded as $RM$ and $RT$, respectively.

Let $\operatorname{ker}(\mathscr{D})$ be the kernel space of operator $\mathscr{D}$. Let $\overline{n}$ denote the congruence of integer $n$ with modulus 3, namely,

(2.4)

$\displaystyle\overline{n}\equiv n(\operatorname{mod}3),~{}~{}~{}~{}~{}\overline{n}=1,2,3.$

For $\bm{u}=(u_{ij})_{3\times 3}$, let $\mathscr{P}_{n}\bm{u}:=(u_{11},u_{22},u_{33})^{\intercal}$ denote the normal stress and $\mathscr{P}_{t}\bm{u}:=(u_{12},u_{13},u_{23})^{\intercal}$ as the shear stress.

Let $\mathcal{T}_{\Box}$ be a cuboid grid of the domain $\Omega\subseteq\mathbb{R}^{3}$. This section constructs finite element subspaces $V_{k,\Box}\subseteq H^{1}(\Omega;\mathbb{R}^{3})$, $\mathcal{U}_{k,\Box}\subseteq H(\operatorname{sym}\operatorname{curl},\Omega;\mathbb{T})$, and $\Sigma_{k,\Box}\subseteq H(\operatorname{div}\operatorname{div},\Omega;\mathbb{S})$ on cuboid grids. It will be proved that these finite element subspaces form an exact discrete divdiv complex of (1.2).

In the endeavor to acquire the exact finite element divdiv complex of (1.2) on $\mathcal{T}_{\Box}$, one of the challenging tasks is to construct the $H(\operatorname{div}\operatorname{div},\Omega;\mathbb{S})$ conforming finite element space $\Sigma_{k,\Box}$. It is not easy to get a unisolvent $H(\operatorname{div}\operatorname{div},\Omega;\mathbb{S})$ conforming finite element space on cuboid grids by strengthening the continuity of the normal component of the divergences of the symmetric matrix valued functions of the $H(\operatorname{div},\Omega;\mathbb{S})$ conforming finite element on cuboid grids as in [12, 6]. To attack this difficulty, the de Rham complex and the divdiv complex which consist of polynomials over hexahedrons are proposed below to determine the shape function space. Besides, the continuity as that of the $H(\operatorname{div}\operatorname{div},\Omega;\mathbb{S})$ conforming finite element space in [12] is imposed sufficiently to construct the DOFs. Then the $H(\operatorname{div}\operatorname{div},\Omega;\mathbb{S})$ conforming finite element space $\Sigma_{k,\Box}$ follows, which is actually a subspace of $H(\operatorname{div}\operatorname{div},\Omega;\mathbb{S})\cap H(\operatorname{div},\Omega;\mathbb{S})$. As a consequence, there is corresponding enhanced smoothness for the remaining finite element spaces $\mathcal{U}_{k,\Box}$ and $V_{k,\Box}$ in the discrete divdiv complex.

For a cuboid partition $\mathcal{T}_{\Box}$, on any cube $K\in\mathcal{T}_{\Box}$ with $K=(x_{a},x_{b})\times(y_{a},y_{b})\times(z_{a},z_{b})$, let

(3.1)

$\displaystyle\begin{split}&b_{K,x}:=\frac{(x-x_{a})(x-x_{b})}{h_{K,x}^{2}},~{}~{}~{}~{}b_{K,y}:=\frac{(y-y_{a})(y-y_{b})}{h_{K,y}^{2}},\\ &b_{K,z}:=\frac{(z-z_{a})(z-z_{b})}{h_{K,z}^{2}}.\end{split}$

Here $h_{K,x}:=x_{b}-x_{a}$, $h_{K,y}:=y_{b}-y_{a}$, and $h_{K,z}:=z_{b}-z_{a}$. Let $b_{K}:=b_{K,x}b_{K,y}b_{K,z}$. And $b_{K,x_{i}}$ is abbreviated as $b_{K,i}$ in the text for $i=1,2,3$. Let $E_{i}(K)\subseteq\mathscr{E}(K)$ be the set of all edges of $K$ which parallel to the $x_{i}$-axis and let $F_{i}(K)\subseteq\mathscr{F}(K)$ be the set of all faces of $K$ that perpendicular to the $x_{i}$-axis. Let $F_{i}(K)\cup F_{j}(K)$ denote the set of all faces of $K$ which are either in $F_{i}(K)$ or in $F_{j}(K)$. Let $F_{i}\subseteq\mathscr{F}$ be the set of all faces of $\mathcal{T}_{\Box}$ that perpendicular to the $x_{i}$-axis. Hereafter, if not special specified, $\{i,j,l\}$ is a permutation of $\{1,2,3\}$.