 # Finite elements for divdiv-conforming symmetric tensors in three dimensions

Two types of finite element spaces on a tetrahedron are constructed for divdiv conforming symmetric tensors in three dimensions. Besides the normal-normal component, another trace involving combination of first order derivatives of stress should be continuous across the face. Due to the rigid of polynomials, the symmetric stress tensor element is continuous at vertices, and on the plane orthogonal to each edge. Hilbert complex and polynomial complexes are presented and several decomposition of polynomial vector and tensors spaces are revealed from the complexes. The constructed divdiv conforming elements are exploited to discretize the mixed formulation of the biharmonic equation. Optimal order and superconvergence error analysis is provided. Hybridization is given for the ease of implementation.

## Authors

##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1. Introduction

Recently we have constructed finite elements for divdiv conforming symmetric tensors in two dimensions . In this paper, we shall continue our study of space

 H(divdiv,Ω;S):={τ∈L2(Ω;S):divdivτ∈L2(Ω)}

and construct corresponding finite element spaces in three dimensions.

It turns out the construction in three dimensions is much harder than that in two dimensions. One obvious reason is the complication of interplay of differential operators ( and ) and matrix operators ( and etc) in three dimensions. The essential difficulty arises from the Hilbert complex

In the divdiv complex in three dimensions, the Sobolev space before

consists of tensor functions rather than vector functions in two dimensions. The subspace of polynomial space with vanishing boundary degree of freedoms is not clear as no finite element spaces for

is known. In two dimensions, we have and finite element spaces for are relatively mature. By analogy, an nonconforming finite element in two dimensions was advanced for discretizing the mixed formulation of the biharmonic equation in [8, 9, 11] in 1960s, while the corresponding nonconforming finite element in three dimensions was constructed to solve a mixed formulation of the linear elasticity until 2011 , rather than the biharmonic equation.

To attack the main difficulty, we present two polynomial complexes and reveal several decomposition of polynomial vector and tensors spaces from the complexes. With further help from the Green’s identity and characterization of the trace, we are able to construct two types of finite element spaces on a tetrahedron. Here we present the BDM-type space below. Let be a tetrahedron and let be a positive integer. The shape function space is simply . The set of edges of is denoted by , the faces by , and the vertices by . For each edge, we chose two normal vectors and . The degrees of freedom are given by

 (1) τ(δ) ∀ δ∈V(K), (2) (n⊺iτnj,q)e ∀ q∈Pk−2(e),e∈E(K),i,j=1,2, (3) (n⊺τn,q)F ∀ q∈Pk−3(F),F∈F(K), (4) (2divF(τn)+∂n(n⊺τn),q)F ∀ q∈Pk−1(F),F∈F(K), (5) (τ,ς)K ∀ ς∈∇2Pk−2(K)⊕sym(x×Pk−2(K;T)), (6) (τn,n×xq)F1 ∀ q∈Pk−2(F1),

where is an arbitrary but fixed face. The degrees of freedom (6) will be regarded as interior degrees of freedom to the tetrahedron , that is the degrees of freedom (6) are double-valued on each face when defining the global finite element space. RT-type space can be obtained by further reducing the index of degree of freedoms by

except the moment with

.

The boundary degree of freedoms (1)-(4) are motivated by the Green’s formulae and the characterization of the trace of . The interior moments of is to determine . Together with , the volume moments can determine the polynomial of degree up to .

We then use the vanished trace and the symmetry of the tesnor. Similarly as the RT and BDM elements , the vanishing normal-normal trace (34) implies the normal-normal part of is zero. To determine the normal-tangential terms, further degrees of freedoms are needed. Due to the symmetry of , it is sufficient to provide additional degrees of freedoms on one face, which are inspired by the RT and BDM elements in two dimensions.

It is arduous to figure out the explicit basis functions dual to the degrees of freedom (1)-(6). Hybridization is thus provided for the ease of implementation. The basis functions of the standard Lagrange element can be used to implement the hybridized mixed finite element methods. The constructed divdiv conforming elements are exploited to discretize the mixed formulation of the biharmonic equation. Optimal order and superconvergence error analysis is provided.

The rest of this paper is organized as follows. In Section 2, we present some operations for vectors and tensors. Two polynomial complexes related to the divdiv complex, and direct sum decompositions of polynomial spaces are shown in Section 3. We derive the Green’s identity and characterize the trace of on polyhedrons in Section 4, and then construct the conforming finite elements for in three dimensions. Mixed finite element methods for the biharmonic equation are developed in Section 5.

## 2. Matrix and vector operations

In this section, we shall survey operations for vectors and tensors. Some of them are standard but some are not unified in the literature. In particular, we shall introduce operators appending to the right side of a matrix for operations applied to columns of a matrix. We will mix the usage of row and column vectors which will be clear in the context.

Given a plane with normal vector , for a vector , we have the orthogonal decomposition

 v=Πnv+ΠFv:=(v⋅n)n+(n×v)×n.

The vector is also on the plane and is a rotation of by counter-clockwise with respect to . Therefore

 ΠFu⋅ΠFv=((n×u)×n)⋅((n×v)×n)=(n×u)⋅(n×v)=Π⊥Fu⋅Π⊥Fv.

We treat Hamilton operator as a column vector and define

 ∇⊥F:=n×∇,∇F:=ΠF∇=(n×∇)×n.

For a scalar function ,

 ∇Fv=ΠF(∇v),∇⊥Fv=n×∇v,

are the surface gradient of and surface , respectively. For a vector function , is the surface divergence and by definition

 ∇F⋅v=∇F⋅(ΠFv)=(n×∇)⋅(n×v)=∇⊥F⋅(n×v).

By the cyclic invariance of the mix product, the surface rot operator is

 ∇⊥F⋅v=(n×∇)⋅v=n⋅(∇×v).

In particular, for , is the plane. Then is the two dimensional operator which is a rotation of the two dimensional divergence operator .

The matrix-vector product can be interpret as the inner product of with the row vectors of . We thus define the dot operator

 b⋅A:=Ab.

Namely the vector inner product is applied row-wise to the matrix. Similarly we can define the cross product row-wise from the left . Here rigorously speaking when a column vector is treat as a row vector, should be used. In most places, however, we will scarify this precision for the ease of notation.

When the vector is on the right of the matrix, the operation is defined column-wise. That is

 A⋅b :=b⊺A=(A⊺b)⊺=(b⋅A⊺)⊺, A×b :=−(b×A⊺)⊺.

By moving the column operation to the right, it is consistent with the transpose operator. For the transpose of product of two objects, we take transpose of each one, switch their order, and add a negative sign if it is the cross product.

For two vectors and matrix , we define the following products

 b⋅A⋅c :=b⊺A⊺c=(b⋅A)⋅c=b⋅(A⋅c), b⋅A×c :=(Ab)×c=(b⋅A)×c=b⋅(A×c), b×A×c :=(b×A)×c=b×(A×c).

By moving the column operation to the right, these operations are associative. That is the ordering of performing the products does not matter.

Apply these matrix-vector operations to the Hamilton operator , we get row-wise differentiation

 ∇⋅A,∇×A,

and column-wise differentiation

 A⋅∇,A×∇.

And we can write the divdiv operator applied to a matrix as

 (7) divdivA=∇⋅A⋅∇.

Denote the space of all matrix by , all symmetric matrix by

, all skew-symmetric

matrix by , and all trace-free matrix by . For any matrix , we can decompose it into symmetric and skew-symmetric part as

 B=sym(B)+skw(B):=12(B+B⊺)+12(B−B⊺).

We can also decompose it into a direct sum of a trace free matrix and a diagonal matrix as

 B=devB+13tr(B)I:=(B−13tr(B)I)+13tr(B)I.

For a vector function , and are standard differential operations. The gradient is a matrix Its symmetric part is defined as

 defu:=sym∇u=12(∇u+(∇u)⊺)=12(u∇⊺+∇u⊺).

In the last identity notation is used to emphasize the symmetry form. Similarly we can define operator for a matrix

 symcurlA=12(∇×A+(∇×A)⊺)=12(∇×A−A⊺×∇).

We define an isomorphism of and the space of skew-symmetric matrices as follows: for a vector

For two vectors , one can easily verify that

 u×v=[u]×v=(mspnu)v.

We will use the following identities which can be verified by direct calculation.

 skw(∇u) =12(mspn∇×u), (8) skw(∇×A) =12mspn[(A⋅∇)⊺−∇(tr(A))], (9) divmspnu =−curlu, (10) curl(uI) =−mspngrad(u).

More identities involving the matrix operation and vector differentiation are summarized in .

## 3. Divdiv complex and polynomial complexes

In this section, we shall consider the divdiv complex and establish two related polynomial complexes. We assume is a bounded and strong Lipschitz domain which is topologically trivial in the sense that it is homeomorphic to a ball. Without loss of generality, we also assume .

### 3.1. The divdiv complex

The complex in three dimensions reads as [1, 13]

 (11) RT\autorightarrow$⊂$H1(Ω;R3)\autorightarrow$devgrad$H(symcurl,Ω;T)\autorightarrow$symcurl$H(divdiv,Ω;S)\autorightarrow$divdiv$L2(Ω)\autorightarrow0.

For completeness, we prove the exactness of the complex (11) following .

###### Theorem 3.1.

Assume is a bounded and topologically trivial strong Lipschitz domain in . Then (11) is a complex and exact sequence.

###### Proof.

Any skew-symmetric can be written as . Assume , it follows from (9) that

 (12) divdivτ=divdivmspnv=−div(curlv)=0.

Since for any smooth skew-symmetric tensor field , we obtain

 divdivH(divdiv,Ω;S)=divdivH(divdiv,Ω;M)=L2(Ω).

For any , it follows from (7) that

 divdivsymcurlτ=12∇⋅(∇×A−A⊺×∇)⋅∇=0.

Hence . For any , it holds from (10) that

We get . As a result, (11) is a complex.

For any , there exists such that

 divσ=curlv=−div(mspnv).

Hence here exists such that

 σ=−mspnv+curlτ.

By the symmetry of , we have Noting that

it follows . Thus

 H(divdiv,Ω;S)∩ker(divdiv)=symcurlH(symcurl,Ω;T).

For any , by the fact that , we have from (8) that

Then

 curl(divτ⊺)=−div(mspndivτ⊺)=−2div(curlτ)=0.

Thus there exists satsifying , which implies

Hence there exists such that . Since is trace-free, we achieve

which means . Therefore the complex (11) is exact. ∎

### 3.2. Polynomial complexes

Given a bounded domain and a non-negative integer , let stand for the set of all polynomials in with the total degree no more than , and denote the tensor or vector version.

###### Lemma 3.2.

The polynomial complex

 (13) RT\autorightarrow$⊂$Pk+2(Ω;R3)\autorightarrow$devgrad$Pk+1(Ω;T)\autorightarrow$symcurl$Pk(Ω;S)\autorightarrow$divdiv$Pk−2(Ω)\autorightarrow0

is exact.

###### Proof.

It follows from (12) that

 divdivPk(Ω;S)=divdivPk(Ω;M)=Pk−2(Ω).

For any , we have . Hence

 −mspn(∇divq)=curl((divq)I)=3curl(∇q)=0,

which means and . We conclude from the fact that

is the identity matrix multiplied by a constant.

For any , there exists satisfying , i.e. . Then

 mspn(∇divv)=−curl((divv)I)=3curl(τ−∇v)=3curlτ∈Pk(Ω;K),

from which we get , and thus . As a result . And we also have

 (15) dimPk(Ω;S)∩ker(divdiv)=16(5k3+36k2+67k+36).

Finally we conclude from (14) and (15). Therefore the complex (13) is exact. ∎

Define operator as

 πRTv:=v(0,0,0)+13(divv)(0,0,0)x.

The following complex is a generalization of the Koszul complex for vector functions.

###### Lemma 3.3.

The polynomial complex

 (16) 0\autorightarrow$⊂$Pk−2(Ω)\autorightarrow$xx⊺$Pk(Ω;S)\autorightarrow$x×$Pk+1(Ω;T)\autorightarrow$x⋅$Pk+2(Ω;R3)\autorightarrow$πRT$RT\autorightarrow0

is exact.

###### Proof.

Since

 tr(x×τ)=2x⊺vspn(skwτ)∀ τ∈L2(Ω;M),

thus (16) is a complex.

For any satisfying , since , there exist and such that . Noting that

 div(τ1x)=x⊺div(τ⊺1)+trτ1=x⊺div(τ⊺1),

we have

 πRT(qx)=πRTv−πRT(τ1x)=0,

which indicates and thus . Hence there exists such that . Taking , we get

 τx=τ1x+xq⊺1x=τ1x+qx=v.

Hence holds.

It holds

 (17) πRTv=v∀ v∈RT.

Namely is a projector. Consequently, the operator is surjective. And we have

 dim(Pk+1(Ω;T)∩ker(x⋅)) =dimPk+1(Ω;T)−dim(Pk+1(Ω;T)x) =dimPk+1(Ω;T)−dimPk+2(Ω;R3)+4 (18) =16(5k3+36k2+67k+36).

For any satisfying , there exists such that . By the symmetry of , it follows

 x×(xv⊺)=x×(vx⊺)⊺=x×τ⊺=0,

which indicates . Then there exists satisfying . Hence . Therefore . And we have

 dim(x×Pk(Ω;S)) =dimPk(Ω;S)−dim(Pk(Ω;S)∩ker(x×)) =dimPk(Ω;S)−dim(Pk−2(Ω)xx⊺) =16(5k3+36k2+67k+36),

which together with (18) implies . Therefore the complex (16) is exact. ∎

Those two complexes (13) and (16) are connected as

 (19)

Unlike the Koszul complex for vectors functions, we do not have the identity property applied to homogenous polynomials. Fortunately decomposition of polynomial spaces using Koszul and differential operators still holds.

Let be the space of homogeneous polynomials of degree . Then by Euler’s formula

 (20) x⋅∇q=kq∀ q∈Hk(Ω).

Due to (20), for any satisfying , we have . And

 (21) Pk(Ω)∩ker(x⋅∇) =P0(Ω), (22) Pk(Ω)∩ker(x⋅∇+ℓ) =0

for any positive integer .

It follows from (17) and the complex (16) that

 Pk+2(Ω;R3)=x⋅Pk+1(Ω;T)⊕RT.

We then move to the space .

###### Lemma 3.4.

We have the decomposition

###### Proof.

Since the dimension of space in the left hand side is the summation of the dimension of the two spaces in the right hand side in (23), we only need to prove that the sum in (23) is the direct sum. For any satisfying and , we have , that is

 (x⋅∇)q=13(divq)x.

By (20), there exist such that . Then

 x⋅∇pi+pi=13divq,

which indicates . Hence it follows from the last identity that , which combined with (21) gives for some constant . Thus , and . This ends the proof. ∎

Finally we present a decomposition of space . Let

 Ck(Ω;S):=symcurlmissingPk+1(Ω;T),C⊕k(Ω;S):=xx⊺Pk−2(Ω).

Their dimensions are

 dimCk(Ω;S)=16(5k3+36k2+67k+36),dimC⊕k(Ω;S)=16(k3−k).
###### Lemma 3.5.

We have

1. for any

2. is a bijection.

###### Proof.

Since and , we get

 (24) divdiv(xx⊺q)=div((x⋅∇+4)qx)=(x⋅∇+3)(x⋅∇+4)q.

Hence property (i) follows from (20). Property (ii) is obtained by writing . Now we prove property (iii). First the dimension of space in the left hand side is the summation of the dimension of the two spaces in the right hand side in (iii). Assume satisfies , which means

 divdiv(xx⊺q)=0.

Thus property (iii) holds from (24) and (22). ∎

For the simplification of the degree of freedoms, we need another decomposition of the symmetric tensor polynomial space which can be derived from the dual complex of polynomial divdiv complexes.

###### Lemma 3.6.

It holds

 (25) Pk(Ω;S)=∇2Pk+2(Ω)⊕sym(x×Pk−1(Ω;T)).
###### Proof.

Noting that the dimension of space in the left hand side is the summation of the dimension of the two spaces in the right hand side in (25), we only need to prove the direct sum.

For any with satisfying , it follows . Applying (21) and (20), we get and . Thus the decomposition (25) holds. ∎

## 4. Divdiv Conforming finite element spaces

We first present a Green’s identity based on which we can characterize the trace of on polyhedrons and give a sufficient continuity condition for a piecewise smooth function to be in . Then we construct the finite element space and prove the unisolvence.

### 4.1. Notation

Let be a regular family of polyhedral meshes of . Our finite element spaces are constructed for tetrahedrons but some results, e.g., traces and Green’s formulae etc, hold for general polyhedrons. For each element , denote by the unit outward normal vector to , which will be abbreviated as for simplicity. Let , , , , and be the union of all faces, interior faces, all edges, interior edges, vertices and interior vertices of the partition , respectively. For any , fix a unit normal vector and two unit tangent vectors and , which will be abbreviated as and without causing any confusions. For any , fix a unit tangent vector and two unit normal vectors and , which will be abbreviated as and without causing any confusions. For being a polyhedron, denote by , and the set of all faces, edges and vertices of , respectively. For any , let be the set of all edges of . And for each , denote by the unit vector being parallel to and outward normal to . Furthermore, set

 Fi(K):=F(K)∩Fih,Ei(F):=E(F)∩Eih.

### 4.2. Green’s identity

We start from the Green’s identity for smooth functions on polyhedrons.

###### Lemma 4.1 (Green’s identity).

Let be a polyhedron, and let and . Then we have

 (divdivτ,v)K =(τ,∇2v)K−∑F∈F(K)∑e∈E(F)(n⊺F,eτn,v)e (26) −∑F∈F(K)[(n⊺τn,∂nv)F−(2divF(τn)+∂n(n⊺τn),v)F].
###### Proof.

We start from the standard integration by parts

 (divdivτ,v)K=−(divτ,∇v)K+∑F∈F(K)(n⊺divτ,v)F=(τ,∇2v)K−∑F∈F(K)(τn,∇v)F+∑F∈F(K)(n⊺divτ,v)F.

We then decompose and apply the Stokes theorem to get

 (τn,∇v)F =(τn,∂nvn+∇Fv)F =(n⊺τn,∂nv)F−(divF(τn),v)F+∑e∈E(F)(n⊺F,eτn,v)e.

Now we rewrite the term

 (n⊺divτ,v)F=(div(τ⋅n),v)F=(divF(τn),v)F+(∂n(n⊺τn),v)F.

Thus the Green’s identity (26) follows by merging all terms. ∎

When the domain is smooth in the sense that is an empty set, the term disappears. When is continuous on edge , this term will define a jump of the tensor.

### 4.3. Traces and continuity across the boundary

We first recall the trace of the space on the boundary of polyhedron (cf. [7, Lemma 3.2] and [17, 15]). Define trace spaces

 H1/2n,0(∂K) :={∂nv|∂K:v∈H2(K)∩H10(K)} ={g∈L2(∂K):g|F∈H1/200(F)∀ F∈F(K)}

with norm

 ∥g∥H1/2n,0(∂K):=infv∈H2(K)∩H10(K)∂nv=g∥v∥2,

and

 H3/2t,0(∂K) :={v|∂K:v∈H2(K),∂nv|∂K=0,v|e=0 for each edge e∈E(K)}

with norm

 ∥g∥H3/2t,0(∂K):=infv∈H2(K)∂nv=0,v=g∥v∥2.

Let for the normal-normal trace, and for the trace involving combination of derivatives.

###### Lemma 4.2.

For any , it holds

 ∥n⊺τn∥H−1/2n(∂K)+∥2divF(τn)+∂n(n⊺τn)∥H−3/2t(∂K)≲∥τ∥H(divdiv).

Conversely, for any and , there exists some such that

 n⊺τn|∂K=gn,2divF(τn)+∂n(n⊺