MATHEMATICS OF COMPUTATION
Volume 80, Number 274, April 2011, Pages 669–695
S 0025-5718(2010)02412-3
Article electronically published on August 26, 2010
DIVERGENCE-FREE FINITE ELEMENTS
ON TETRAHEDRAL GRIDS FOR k ≥ 6
SHANGYOU ZHANG
Abstract. It was shown two decades ago that thePk-Pk−1 mixed element
on triangular grids, approximating the velocity by the continuousPk piecewise
polynomials and the pressure by the discontinuousPk−1 piecewise polynomi-
als, is stable for all k ≥ 4, provided the grids are free of a nearly-singular
vertex. The problem with the method in 3D was posted then and remains
open. The problem is solved partially in this work. It is shown that thePk-
Pk−1 element is stable and of optimal order in approximation, on a family of
uniform tetrahedral grids, for allk ≥ 6. The analysis is to be generalized to
non-uniform grids, when we can deal with the complicity of 3D geometry.
For the divergence-free elements, the ﬁnite element spaces for the pressure
can be avoided in computation, if a classic iterated penalty method is applied.
The ﬁnite element solutions for the pressure are computed as byproducts from
the iterate solutions for the velocity. Numerical tests are provided.
1. Introduction
Rewriting the Navier-Stokes or the Stokes equations in the weak variational
forms, the primitive unknowns, the velocity and the pressure, belong to Sobolev
spacesH1 andL2, respectively. Naturally, a ﬁnite element method would bethePk-
Pk−1 element which approximates the velocity in anH1-subspace of continuousPk
piecewise polynomials (C0-Pk) and approximates the pressure in anL2-subspace of
discontinuous Pk−1 piecewise polynomials (C−1-Pk−1). This is a truly conforming
element as the incompressibility condition is satisﬁed pointwise and the discrete
solution for the velocity is a projection within the space of divergence-free functions.
A fundamental study on the method was done by Scott and Vogelius ([11, 12]) that
the method is stable and consequently of the optimal order on 2D triangular grids
for anyk ≥ 4, provided that the grids have no nearly-singular vertex. A 2D vertex
of a triangulation is singular if all edges meeting at the vertex form two cross lines;
see Figure 1. For k ≤ 3, Scott and Vogelius showed that thePk-Pk−1 element
would not be stable, and may not produce approximating solutions on general 2D
triangular grids in [11, 12]. What is this magic numberk in 3D? Scott and Vogelius
posted this question explicitly after discovering thatk = 4 in 2D. The problem has
remained open for more than 20 years.
The geometry of the 3D tetrahedral grids is much more complicated than that of
2D. By adding or moving a few edges and vertices locally, one can easily eliminate
Received by the editor June 18, 2008 and, in revised form, January 25, 2010.
2010 Mathematics Subject Classiﬁcation.Primary 65N30, 76M10, 76D07.
Key words and phrases. Mixed ﬁnite elements, Stokes equations, divergence-free element,
tetrahedral grids.
c© 2010 American Mathematical Society
Reverts to public domain 28 years from publication
669

670 SHANGYOU ZHANG
 
 
 
 
  
@
@
@
@
@@
 
 
 
 
  
@
@
@
@
@@
A







B




S
S
S
S
S
S
SS
QQQQQQ







C
Figure 1. Singular vertices (A and B) and a nearly-singular ver-
tex (C,w h e nC → A), in 2D.
singular vertices in 2D; see Figure 1. When a triangulation is singular-vertex free, it
is shown by Scott and Vogelius [11, 12] that the divergence of aC0-Pk vector space
is exactly the space ofC−1-Pk−1 modulus a constant. Following this approach, we
found previously that this is true on Hsieh-Clough-Tocher tetrahedral grids ([17])
for allk ≥ 3 in 3D also. For general tetrahedral grids, it is challenging to identify
all the singular vertices and edges. For example, when doing multigrid reﬁnements
on tetrahedral grids (cf. [16]), a known type of singular edges (all face triangles
meeting at the edge fall into two planes) and singular vertices (all face triangles
meeting at the vertex fall into three planes) cannot be avoided. To extend the
Scott-Vogelius result to 3D while avoiding the technical details on the geometry
of singular vertices, we limit this research on a family of uniform grids, shown in
Figure 2. We will show that thePk-Pk−1 element is stable and provides the optimal
order solutions, for allk ≥ 6. When a classic iterated penalty method ([6, 3, 4, 14])
is used here, we only need to solve a vector-Laplacian equation for the velocity with
an iteration number independent of grid size. In such a case, the mixed element
is reduced to a single element, and the pressure is computed as a byproduct. This
research is still far away from answering the question on the magic numberk in 3D
proposed by Scott and Vogelius. Since we limit our work on the uniform grids, the
magic k may be greater than 6. As we requirek ≥ 6 in our constructional proof,
the magick could be less than 6 as well, though unlikely; see Corollary 3.1 and the
numerical result following that. We note that for the continuous pressure version
of the Pk-Pk−1 element (k ≥ 2) on tetrahedral grids, the analysis is done in [2],
extending the Taylor-Hood element [10].
The rest of the paper is organized as follows. In Section 2, we deﬁne thePk-
Pk−1 element. In Section 3, we will prove the stability of thePk-Pk−1 element on a
uniform grid, and show the optimal order of convergence. In Section 4, we provide
some numerical results.
2. The Pk-Pk−1 element
In this section, we shall deﬁne the Pk-Pk−1 ﬁnite element for the stationary
Stokes equations. The resulting linear systems are guaranteed to have a unique
solution, i.e. the (reduced) inf-sup condition always holds for such a divergence-
free ﬁnite element pair. The classic iterated penalty method ([6, 3, 4, 14]) can be
applied where the mixed element is reduced to a single divergence-free element.

DIVERGENCE-FREE FINITE ELEMENTS ON TETRAHEDRAL GRIDS 671
Figure 2. The ﬁrst three levels (n =1 ,2,4) of grids, Ωh.
Wesolve amodelstationaryStokesproblem: Find functions u(theﬂuid velocity)
and p (the pressure) on a domain of unit cube Ω = (0,1)3 such that
(2.1)
−Δu+∇p = f in Ω,
divu = 0 in Ω ,
u = 0 on ∂Ω,
where f is the body force. The standard variational form is: Findu ∈ H1
0(Ω)3 and
p ∈ L2
0(Ω) such that
(2.2)
a(u,v)+ b(v,p)=( f,v) ∀v ∈ H1
0(Ω)3,
b(u,q)=0 ∀q ∈ L2
0(Ω).
Here H1
0(Ω)3 is the Sobolev space (cf. [5]) with zero boundary trace,L2
0(Ω) is the
L2 space with zero mean value, i.e.,L2(Ω)/R = {p ∈ L2 |
∫
Ω p =0 },a n d
a(u,v)=
∫
Ω
∇u·∇v dx,
b(v,p)= −
∫
Ω
divu pd x,
(f,v)=
∫
Ω
fv dx.
Let Ωh be a family of uniform tetrahedral grids on Ω depicted in Figure 2:
Ωh = {K | K is a tetrahedron with size|K|≤ h}.
Then we deﬁne thePk-Pk−1 mixed element spaces by
Vh,k =
{
uh ∈ C(Ω) | uh|K ∈ Pk(K)3 ∀K ∈ Ωh and uh|∂Ω =0
}
⊂ H1
0(Ω)3,
(2.3)
Ph = {divuh | uh ∈ Vh,k}⊂ L2
0(Ω).(2.4)
It is widely known that the pointwise divergence-free mixed method is too compli-
cated and not practical; cf. [4]. Very little work has been done on this method; cf.
[1, 7, 8, 9, 11, 17, 18]. The resulting system of ﬁnite element equations for (2.2) is:
Find uh ∈ Vh,k and ph ∈ Ph such that
(2.5)
a(uh,v)+ b(v,ph)=( f,v) ∀v ∈ Vh,k,
b(uh,q)=0 ∀q ∈ Ph.

672 SHANGYOU ZHANG
The linear system of equations (2.5) always has a unique solution, in the divergence-
free element method; cf. [18]. We note that Ph in (2.4) is a proper subspace of
traditional C−1-Pk−1 ﬁnite element space. We will characterize it in detail below.
Letting q =d i vuh in (2.5), we still have the (pointwise) divergence-free property
for the ﬁnite element solution
(2.6)
∫
Ω
(divuh)2dx = b(uh,q)=0 .
By (2.6), the unique solutionuh of (2.5) is divergence-free ([10, 4, 3, 18]). It is, in
fact, thea(·,·) orthogonal projection from the divergence-free spaceZ to a subspace
Zh, deﬁned by,
Z :=
{
v ∈ H1
0(Ω)3 | divv =0
}
,(2.7)
Zh := {v ∈ Vh,k | divv =0 }.(2.8)
As Ph may be a proper subspace of discontinuous, piecewise polynomials of degree
(k−1) or less, it may be diﬃcult to ﬁnd a nodal basis forPh in some cases. But on
the other side, it is the special interest of the divergence-free element method that
the spacePh can be omitted in computation and the discrete solutions approximat-
ing the pressure function in the Stokes equations can be obtained as byproducts, via
the iterated penalty method. This does not only simplify the coding work, but also
it avoids the diﬃculty of solvingnon-positive deﬁnite systems of linear equations,
encountered in typical mixed element methods. We refer to [6, 4, 3, 14, 18] for the
iterated penalty method.
3. Stability and convergence
In this section, we will prove the inf-sup condition (3.69), i.e., the stability of
the divergence-freePk-Pk−1 mixed element. The analysis is done by construction,
based on the unit cube domain Ω and the uniform grids Ωh, except Lemma 3.1.
The convergence follows the stability routinely.
Lemma 3.1.For anyq ∈ Ph (deﬁned in(2.4)), k ≥ 3, there is a functionv1 ∈ Vh,3
(deﬁned in(2.3)) such that
(3.1)
∫
K
divv1 =
∫
K
q ∀K ∈ Ωh, and ‖v1‖H1(Ω)3 ≤ C‖q‖L2(Ω).
Proof. For any q ∈ Ph, by the inf-sup condition for the continuous functions
(cf. [10]) there is auq ∈ H1
0(Ω)3 such that
divuq(x,y,z )= q(x,y) a.e. for ( x,y,z ) ∈ Ω
and
‖uq‖H1 ≤ C‖q‖L2.
We modify the Lagrange interpolation operator slightly to deﬁne a “Fortin op-
erator” (see [4]):
Ih : C(Ω)∩H1
0(Ω)3 → Vh,3, Ih : uq ↦→ Ihuq,
Ihuq(ai)= uq(ai) at all nodes except the four internal face nodes,
∫
(∂K)i
Ihuqdx =
∫
(∂K)i
uqdx,i =1 ,2,3,4,

DIVERGENCE-FREE FINITE ELEMENTS ON TETRAHEDRAL GRIDS 673









J
J
J
J
J
J
J
JJ

@
@@
A
A
A
A
AA
ai :
s
s
s
s
sssss
s
s
ss
s
s
s
s
s
s
ss
s
s
s









J
J
J
J
J
J
J
JJ

@
@@
A
A
A
A
AA
q
q
q
q
qqqqq
q
q
qq
q
q
q
q
q
q
qq
q
q
q
cc
c
bi :
Figure 3. P3 Lagrange nodes: bi. Removing four inner-face nodes.
where Ihuq(bi) (see Figure 3) is chosen so that the integral on each of the four
face triangles matches that ofuq. We note that an averaging interpolation can be
adopted if the functionuq is not continuous, as usual; see [13]. Also follow, for
example, [13], it is standard to show the stability of such an interpolation operator
by scaling:
‖Ihuq‖H1 ≤ C‖uq‖H1.
The interpolant also preserves the divergence elementwise:∫
K
divv1dx =
∫
∂K
v1 ·ndx =
∫
∂K
uq ·ndx =
∫
K
divuqdx =
∫
K
qdx.
We note that the above analysis in deﬁningIhuq ∈ Vh,3 is well known in showing
the stability ofP3-P0 element in 3D; cf. [17]. □
After matching the integral values ofq elementwise by divv1,w en e x tm a t c ht h e
vertex-values ofq −divv1.
Lemma 3.2. For any q ∈ Ph deﬁned in (2.4) such that
∫
K q =0 ∀K ∈ Ωh,
k ≥ 3, there is a functionv2 ∈ Vh,3 such that
divv2(aK
i )= q(aK
i ) ∀K ∈ Ωh,(3.2)
∫
K
divv2 =0 ∀K ∈ Ωh,(3.3)
‖v2‖H1(Ω)3 ≤ C‖q‖L2(Ω).(3.4)
Here aK
i , 1 ≤ i ≤ 4, are the four vertices of elementK.
Proof. Let q =d i vwh for some wh ∈ Vh,k, k ≥ 3. From Figure 4, there are six
types of vertices in Ωh:
Type (a): Corner vertices shared by 2 tetrahedra,B,C,D,E,F,H in Figure 4,
Type (b): Corner vertices shared by 6 tetrahedra,A and G in Figure 4,
Type (c): Mid-edge vertices shared by 4 tetrahedra,I, P and R in Figure 4,
Type (d): Mid-edge vertices shared by 8 tetrahedra,J in Figure 4,
Type (e): Mid-face vertices shared by 12 tetrahedra,L, N and Q in Figure 4,
Type (f): Internal vertices shared by 24 tetrahedra,M in Figure 4,
For a Type (a) boundary vertex, such asB in Figure 4, the vector ﬁeldwh
vanishes on the four boundary faces meeting atB; it follows that in each of the
two tetrahedra sharing the vertexB, wh vanishes along the three edges meeting at

674 SHANGYOU ZHANG
E
F
B
A
G
H
D
C
E
F
B
A
G
H
D
C
I
J
M
P
N
S
R
Q
L
Figure 4. The interior and boundary vertices of Ωh.
B (in BASGF , for instance,wh vanishes alongBA, BF,a n dBG). This implies
that divwh =0a t B, so that we do not need any construction ofv2 in order to
meet the requirement (3.2), becauseq|GAFB (B)= q|GACB(B)=0 .
For a Type (c) vertex, similarly, all four tetrahedra meeting at the vertex have
three boundary edges. Therefore, q(aKj
i )=d i vwh(aKj
i ) = 0 at 4 tetrahedraKj,
sharing such a boundary vertex.
For a Type (b) vertex such asA in Figure 4, there are six tetrahedra {Kj}
sharing the vertex. We deﬁne a vector functionv2,(b) ∈ P3
3 ∩ C0(∪Kj) such that
v2,(b) = 0 at all Lagrange nodes except nodes on the diagonal edge of the cube
formed by the six tetrahedra, i.e., nodesa and b in Figure 5. As v2,(b) has three
components, there are in total 6 degrees of freedom for such av2,(b),a tn o d e sa
and b. Note that, as wh|∂Ω = 0, the gradient ofwh at A are the same on two
tetrahedra sharing a ﬂat boundary. That is,
q|AGEH(A)=q|AGHD(A),q |AGDC(A)=q|AGCB(A),q |AGBF (A)=q|AGFE (A).
The three values ofq(A) at a boundary vertexA would be matched by three degrees
of freedom ofv2,(b) while the other three degrees of freedom ofv2,(b) would make
∇v2,(b) = 0 at the opposite vertexG.
Let us give an explicit construction of v2,(b). Without loss of generality, let
ABCDEFGH be the unit cube at the origin. Let a “derivative nodal” basis func-
tion φ(b)(x,y,z )a tA be the continuous piecewiseP3 function which has nodal value
0 at all Lagrange nodes except two diagonal nodesa and b (see Figure 5), so that
∇φ(b)(G)=
⎛
⎝
0
0
0
⎞
⎠, but ∇φ(b)(A)=
⎧
⎪⎪⎪⎪⎨
⎪⎪⎪⎪⎩
(
100
)T
on AGEH ∪AGHD,
(
001
)T
on AGDC ∪AGCB,
(
010
)T
on AGBF ∪AGFE.

DIVERGENCE-FREE FINITE ELEMENTS ON TETRAHEDRAL GRIDS 675
A
B
D
C
E
F
H
G
A
E
GF
a
b
c
A
E
G
H
a
b
Figure 5. The vertices and the Lagrange interpolation nodes.
We can deﬁne (not unique) this continuousP3 nodal basis by (see Figure 4),
(3.5) φ(b)(x,y,z )=
⎧
⎪⎪⎪⎪⎪⎪⎪⎪⎨
⎪⎪⎪⎪⎪⎪⎪⎪⎩
x(1−z)(1+ y −2z)o n AGEH,
x(1−y)(1+ z −2y)o n AGHD,
z(1−y)(1+ x−2y)o n AGDC,
z(1−x)(1+ y −2x)o n AGCB,
y(1−x)(1+ z −2x)o n AGBF,
y(1−z)(1+ x−2z)o n AGFE.
The functionv2,(b) to be constructed is
(3.6) v2,(b) = φ(b)(x,y,z )
⎛
⎝
q|AGEH(A)
q|AGBF (A)
q|AGDC(A)
⎞
⎠.
By the construction divv2,(b) has zero nodal values at all vertices of 6Kj, except
at vertexA, where the six values match that ofq. By the equivalence of norms on
the unit cube for piecewise polynomials, we have
(3.7) |q(A)|≤ Ch−3/2‖q‖L2(⋃Kj).
On the otherside, we used scaled derivatives todeﬁnev2,(b) and we get thefollowing
bound:
|v2,(b)|H1(⋃Kj)3 ≤ C|divv2,(b)|L2(⋃Kj) ≤ Ch3/2|divˆv2,(b)|L2((0,1)3)
≤ Ch3/2|q(A)|≤ C‖q‖L2(⋃Kj).(3.8)
We note that due to the uniform grid, we can compute the constants in (3.7) and
(3.8). For example, by (3.5) and (3.6), we can obtain
|v2,(b)|H1(AGEH)3 =
⏐⏐q|AGEH(0,0,0)
⏐⏐√
3|x(1−z)(1+ y −2z)|H1(AGEH)
=
⏐⏐q|AGEH(0,0,0)
⏐⏐1√
35
=2
√
3‖divv2,(b)‖L2(AGEH).
But the constants in (3.7) and (3.8) would depend on the polynomial degreek.
The constructedv2,(b) satisﬁes (3.2) and (3.4), but not (3.3), i.e.,
∫
Kj
divv2,(b) ⁄=
0. By the divergence theorem, the integral on the whole cube formed by the 6
tetrahedra is zero. After correcting the integrals on 5 of the 6 tetrahedron by

676 SHANGYOU ZHANG
functions supported inside two tetrahedra each time, the last integral on the sixth
integral would be zero also. First, we deﬁnev2,(b1) by
(3.9) v2,(b1) =
{
c0nAGEφAGEH on tetrahedronAGEH,
d0nAGEφAGEF on tetrahedronAGEF,
where nAGE is the outward normal to the faceAGE on tetrahedron AGEH and
φAGEH is aP3 polynomial identically zero on the three faces ofAGEH except face
triangle AGE. φAGEH is zero on all Lagrange nodes exceptc; cf. Figure 5. For
example, ifAGEH is the unit cube as in (3.6),φAGEH = x(1−z)(z−y). In (3.9),
c0 is chosen so that
(3.10)
∫
AGEH
divv2,(b1)dx =
∫
AGEH
divv2,(b)dx.
In (3.9)d0 is chosen so thatv2,(b1) is continuous on the interface. We note thatc0
can always be found to satisfy (3.10) asnAGE ·∇φAGEH is strictly positive inside
AGEH for the third degree polynomialφAGEH. By a scaling argument, we have
also that
(3.11) ‖v2,(b1)‖H1(Ω)3 ≤ C‖v2,(b)‖H1(Ω)3 ≤ C‖q‖L2(Ω).
Next, we repeat the process on the two tetrahedraAGFE and AGFB to deﬁne
v2,(b2) so that
∫
AGEH
divv2,(b2)dx =
∫
AGEH
(
divv2,(b) −divv2,(b1)
)
dx.
It follows by the construction that
(3.12) ‖v2,(b2)‖H1(Ω)3 ≤ C‖v2,(b1)‖H1(Ω)3 +C‖v2,(b)‖H1(Ω)3 ≤ C‖q‖L2(Ω).
Repeatedly, we obtainv2,(bi), i =1 ,2,..., 5. We note that after we deﬁnev2,(b5),
the integral of the divergence of the diﬀerencev2,(b) −∑v2,(bi) h a st ob ez e r oo n
the sixth tetrahedron as it is zero on the seventh tetrahedron which is also the ﬁrst
tetrahedron. Let ˜v2,(b) = v2,(b) − ∑v2,(bi).T h e n˜v2,(b) satisﬁes (3.3) and (3.4),
and its divergence matchesq at 6 vertices of 6 tetrahedra atA while being zero at
all other vertices. By symmetry, we can construct such a˜v2,(b) at the other vertex
G.
I
L
U G
J
H
Q
M
T S
N
P
E
F
B
A
G
H
D
C
I
J
M
P
N
S
R
Q
L
R′ D′
Figure 6. 8 tetrahedra meeting at a vertexJ.

DIVERGENCE-FREE FINITE ELEMENTS ON TETRAHEDRAL GRIDS 677
For a Type (d) vertex, for example,J in Figure 4, we construct a˜v2,(d).T h i s
time, we have two internal edges meeting atJ, on which we have 6 degrees of
freedom to match divv2,(d) with q at the 8 vertices of 8 tetrahedra meeting atJ.
Similar to (3.5) and (3.7), we deﬁne one part ofv2,(d) as follows (cf. Figure 4).
divv2,(d1)|ILJ(J)= q|ILJ(J)= q|IHJ (J),(3.13)
divv2,(d1)|HPJ (J)= q|NPJ (J)= q|IHJ (J),(3.14)
divv2,(d1)|MQLJ (J)= q|MJLQ (J) ⁄= q|MJQN (J).(3.15)
To do so, we repeat the construction ofv2,(b). Next, on an internal edgeJM,w e
deﬁne a nodal basis functionφ(d) like (3.5):
φ(d) =
⎧
⎪⎪⎪⎨
⎪⎪⎪⎩
(1−y)(z −x)(−y)o n MJNG,
(y −x)(1−z)(−y)o n MJGL,
(1+ x−z)(y)(−y)o n MJLQ,
(1+ x−y)(z)(−y)o n MJQN,
assuming thatM is the origin andMN is an edge in they direction of length 1.
We construct the second part ofv2,(d) by
v2,(d2) = φ(d)
⎛
⎝
q|MJNG (J)−(q|MJLQ (J)−q|MJQN (J))
q|MJNG (J)
q|MJNG (J)−(q|MJLQ (J)−q|MJQN (J))
⎞
⎠.
Then we letv2,(d) = v2,(d1) +v2,(d2).d i vv2,(d) matches q at 7 vertices atJ, except
div(v2,(d1) +v2,(d2))|MJGL (J)= q|MJNG (J)−q|MJLQ (J)+ q|MJQN (J).
Will divv2,(d)|MJGL (J)= q|MJGL (J)? The answer is yes. As continuousP3 func-
tions, the gradients of three componentswh at J are
gradwh,i =
⎧
⎪⎪⎪⎪⎪⎪⎪⎨
⎪⎪⎪⎪⎪⎪⎪⎩
(
0 wi,1 0
)T
on MJNG,
(
00 wi,1
)T
on MJGL,
(
wi,2 0 wi,1
)T
on MJLQ,
(
wi,2 wi,1 0
)T
on MJQN,
where wi,j are constants. Then
divwh(J)=
⎧
⎪⎪⎪⎨
⎪⎪⎪⎩
w2,1 on MJNG,
w3,1 on MJGL,
w1,2 +w3,1 on MJLQ,
w1,2 +w2,1 on MJQN,
i.e.,
divwh|MJGL (J)+div wh|MJQN (J)=d i vwh|MJNG (J)+div wh|MJLQ (J).
Hence, as divv2,(d) matches the 7 values ofq =d i vwh at J, it matches the eighth
value divwh|MJGL (J). Finally, we correct the perturbation of v2,(d) on the 8
tetrahedra by 7 bubble functionsv2,(dj) to obtain a ˜v2,(d) to preserve condition
(3.3), as we did for˜v2,(b) by {v2,(bi)}1≤i≤5.

678 SHANGYOU ZHANG
For a Type (e) vertex, say, the mid-face vertexN in Figure 4. The additional
directional derivative ofwh at N of internal edgeNM will give us aP3 nodal basis
(cf. Figure 6)
φ(e) =
⎧
⎪⎪⎪⎪⎪⎪⎪⎪⎨
⎪⎪⎪⎪⎪⎪⎪⎪⎩
(1−y)(y −x)(y)o n MNGS,
(1−y)(y −z)(y)o n MNJG,
(1+ x−y)(y −z)(y)o n MNQJ,
(1+ x−y)(y)(y)o n MND′Q,
(1+ z −y)(y)(y)o n MNR′D′,
(1+ z −y)(y −x)(y)o n MNSR ′,
assuming M is the origin andMN is a unit edge in they direction. Let v2,(e1) be
φ(e)
(
−q|MNGS (N)
)
⎛
⎝
0
1
0
⎞
⎠.
Then, as wh vanishes on boundary triangles GNJ and GNS,w eh a v et h a t
divwh|MNGS (N)=d i vwh|MNJG (N), and that divv2,(e1) matches q at N on the
two tetrahedra. Next, viewing the bottom two cubes (having face squaresCSNR
and NRDP , respectively) belowN together, N is a typeJ node. So, by the con-
struction ofv2,(d), we match nodal values of (q−divv2,(e1))a tN by a 5-dimensional
space to get av2,(e2). Again, viewing the two cubes (having face squaresRDPN
and PNJH , respectively) behind N together, this also makesN at y p eJ node
(a vertical mid-edge Type (d) node). We can deﬁne anotherv2,(e3) to match its
divergence with (q −divv2,(e1) −divv2,(e2))a tN. Therefore, the divergence of
v3,(3) = v2,(e1) +v2,(e2) +v2,(e3)
matches q at N, on all 12 tetrahedra. Unlike earlier cases, we do have (enough
count) 12 degrees of freedom atN in {vh} (3 components and 4 internal edges),
but {(divvh)(N)} is only of dimension 8. On the other side, the twelve values ofq
at N would also form an eight-dimensional vector space, becauseq =d i vwh and
we have the following four constraints:
q|MNGS (N)= q|MNGL (N),(3.16)
q|D′NDR(N)= q|D′NDP (N),(3.17)
q|R′NMS (N)−q|R′NSR(N)= −q|R′NRQ(N)+ q|R′ND′M(N),(3.18)
−q|QNJM (N)+ q|QNMD ′(N)= q|QND′P(N)−q|QNPJ (N).(3.19)
Repeating the process (3.5)–(3.12), after correcting the integral of divergence of
divv2,(e) on 12 tetrahedra by 11 bubble functions, we would obtain a˜v2,(e) for the
lemma.
For a Type (f) node, at an internal vertexM in Figure 4, we have 14 internal
edges and 24 tetrahedra connected to the vertex. These 3 × 14 = 42 degrees
of freedoms for v2,(f) will make the divergence of it matchq values at M on 24
tetrahedra. Here {q|Ki(N)} is a dimension 18 vector space, not of 24 dimensions.
There are 6 constraints similar to (3.18) and (3.19), around the 6 square-diagonal
edges meeting at N; cf. Figure 6. To construct v2,(f), we ﬁrst view M as two
overlapping Type (e) vertices with 4 squares on the left and 4 on the right toM.
Then we separate the eight squares meetingM into two groups, 4 on top and 4 at

DIVERGENCE-FREE FINITE ELEMENTS ON TETRAHEDRAL GRIDS 679
the bottom, in order to use the construction for a type (e) boundary vertex. Of
course, we can constructv2,(f) directly by giving its explicit deﬁnition as we did
for v2,(d) and v2,(e). Repeatedly, we correctv2,(e) to get ˜v2,(f) to preserve (3.3).
The lemma is proved by letting
v2 =
∑
2 Type (b) vertices
˜v2,(b) +
∑
4(n−1) Type (d) vertices
˜v2,(d)
+
∑
6(n−1)2 Type (e) vertices
˜v2,(e) +
∑
(n−1)3 Type (f) vertices
˜v2,(f),
where n is the number of cubes in one direction. □
After we match the element integrals and the vertex values ofq ∈ Ph, we will
next matchq pointwise on each edge within each element.
Lemma 3.3. For anyq ∈ Ph deﬁned in(2.4) such that
∫
K q =0 ∀K ∈ Ωh, k ≥ 6
and q vanishes at all vertices of gridΩh, there is a functionv3 ∈ Vh,k such that
divv3|EK
i
= q|EK
i
∀K ∈ Ωh,(3.20)
∫
K
divv3 =0 ∀K ∈ Ωh,(3.21)
‖v3‖H1(Ω)3 ≤ C‖q‖L2(Ω).(3.22)
Here EK
i , 1 ≤ i ≤ 6, are the six edges of tetrahedronK.
Proof. Let wh ∈ Vh,k such that divwh = q for aq satisfying the lemma conditions.
We will construct av3 matching its divergence with divwh at all edges. We start
with an edgeEA of triangleEAG; see Figure 8. As in Lemma 3.2, we ﬁrst construct
a v3,1 in Vh,k for wh ∈ Vh,k, for all k ≥ 4, matching q at edge EA. Then, we
correct the integral of divv3,1 by av3,0 supported on two tetrahedraFAGE and
EAGH in Figure 8. Regardless of the polynomial degree k in the lemma, the
polynomial degree for v3,0 can be chosen exactly 6 for allk ≥ 4. The reason is
that in order to correct the elementwise divergence-free condition (3.21) while not
perturbing the divergence on the edges shared by two neighboring tetrahedra, as
we did in (3.10)–(3.12), we need degree 6 “bubble” polynomials which have an
internal-face degree of freedom, shown in Figure 9.
r
r
r
r
r
r
r
r
rrrrd
d
d
d
d
d
d
d
dddd
@
@
@
@
@
@
@
@
@
rrrrrrr r
r
r
r
r
r
r
r
r
r
r
r
r
r
r
9“ rd” Lagrange nodes to
9“ 6” Hermit nodes.
6 6 6 6
-
-
-
  	
  	
@
@
@
@
@
@
@
@
@
rrrrrrr r
r
r
r
r
r
r
r
r
r
r
r
r
r
r
Figure 7. Change {λi =1 ,λj > 0} Lagrange nodes to Hermit
nodes, forP6 elements.

680 SHANGYOU ZHANG
A
B
D
C
E
F
H
G
A
E
GF
A
E
G
H
Figure 8. Modiﬁed Lagrange interpolation nodes forP7.
To match the divergence ofwh at edgeAE (see Figures 8 and 9) we replace some
of the standard Lagrange interpolation nodes on the face triangleAEG by some
edge-normal derivatives on the face triangle. Here we replace one loop of Lagrange
nodes of the standardPk element on one face triangle by Hermit nodes on the three
edges of the triangle, shown in Figure 7 and Figure 8. We show next that thePk
element is well deﬁned this way, fork ≥ 4. Let v ∈ Pk deﬁned on tetrahedron
AEGH so that all interpolation values are zero. Let the restriction ofv on triangle
AEG be vk.L e tLAE =0 ,LEG =0a n dLGA = 0 be the equations for three lines
AE, EG and GA, respectively. Sincevk has (k +1) zero points on the three lines,
we have
(3.23) vk = LAELEGLGAvk−3, for somevk−3 ∈ Pk−3(AEG)2.
Let nLAE be the unit normal vector toAE inside plane EAG.A s ∂vk/∂nLAE
has (k −2) zero points on the lineAE (see Figures 8 and 9)vk−3|AE ≡ 0.
vk = L2
AELEGLGAvk−4, for somevk−4 ∈ Pk−4(AEG)2.
Again, as∂vk/∂nLEG has (k −3) zero points on the lineEG,a n d∂vk/∂nLGA has
(k −4) zero points on the lineGA, it follows that
vk = L2
AEL2
EGL2
GAvk−6, for somevk−6 ∈ Pk−6(AEG)2.
Finally, asv =0a t( k −4)(k −5)/2 Lagrange nodes interior to triangleAEG,w e
conclude thatvk−6 ≡ 0a n dv|AEG ≡ 0. Therefore,
(3.24) v = LAEGwk−1, for somewk−1 ∈ Pk−1(AEGH)3.
Here LAEG = 0 is an equation for the planeAEG. As the rest of the Lagrange
interpolationpointsarenotaltered, wk−1 =0a t(k+2)(k+1)k/6standardLagrange
nodes forPk−1 in 3D (see Figure 8) we conclude thatwk−1 =0a n dv =0 .
Now we are ready to prove the lemma. For each internal triangle, exactly one
face triangle of each of two tetrahedra sharing this internal triangle are on the same
plane, due to special structure of the uniform grid. For example, the edgeEG of
internal triangleEGA is on the planeEFGH of two face trianglesEFG and EGH
of the two tetrahedraEFGA and EGHA; cf. Figure 4. For internal triangleQNP,
thetwosharing tetrahedrahaveedge QN on theirtwo facetriangles plane,IQRNJ .
If such an edge is on the boundary, then we lose all internal edge degrees of freedom

DIVERGENCE-FREE FINITE ELEMENTS ON TETRAHEDRAL GRIDS 681
Figure 9. Modiﬁed Lagrange interpolation nodes forP4, P5, P6
and P7.
for divwh on one side of the edge, for example, divwh|AFE (x)=d i vwh|AFB (x)
for x on edgeAF; see Figure 4.
Let us ﬁrst try to construct aP4 polynomial v3 at two tetrahedraGFEA and
GHEA sharing an internal triangleAEG (cf. Figure 8) to match divv3 and divwh
on the three edges of triangleAEG. Let us try to deﬁne aP4 “derivative nodal
basis” shown in Figure 9, which is 0 on the 6 outside face triangles ofGFEA
and GHEA and has a normal derivative 1 inside the faceEAG and two normal
derivatives 0 on another edge. For example, onEAGH,w eh a v e
φ3,1(x,y,z )= x(1−z)(z −y)(c1 +c2x+c3y +c4z)
with 4 constants. But one of them is determined by the internalP4 Lagrange node,
inside the tetrahedron. The other three constants would be determined by three
normal derivatives, 2 on one edge, 1 on another, 0 on the third edge, shown in
Figure 9. We have three choices, EA, EG or AG, for the two-derivative edge,
where the 2 normal derivatives inside triangleEAG are 0. This gives us three such
“nodal basis” functions:
(3.25)
φ3,1 =
{
27
4 x(1−z)(z −y)x on EAGH,
27
4 y(1−z)(z −x)x on EAGF,
φ3,2 =
{ 27√
2x(1−z)(z −y)(1−z)o n EAGH,
27√
2y(1−z)(z −x)(1−z)o n EAGF,
φ3,3 =
{ 27
2
√
2x(1−z)(z −y)(z −x)o n EAGH,
27√
2y(1−z)(z −x)(z −x)o n EAGF.
Since a vh function has three components, with 3φ3,i we can have a 3× 3=9
dimensional subspace
(3.26) {vh} =s p a n
⎧
⎨
⎩φ3,iej | e1 =
⎛
⎝
1
0
0
⎞
⎠, e2 =
⎛
⎝
0
1
0
⎞
⎠, e3 =
⎛
⎝
0
0
1
⎞
⎠
⎫
⎬
⎭;
see (3.27) below. However, we need a dimension 10{divvh} for the 10 degrees
of internal-edge freedom ofq on the two sides of triangleEAG. In fact, we need
a dimension 12{divvh} subspace for a general grid. But we have a special grid
here that every triangle has precisely one singular edge, where the two neighboring
tetrahedra have one common face plane. In this case, the edgeEG is a singular

682 SHANGYOU ZHANG
edge asEGFA and HEGA each have a triangle on planez = 1. When a singular
edge is on the boundary, q is continuous on it. In particular, when k =4 ,w e
have q|EAGF (1
3, 1
3,1) = q|EAGH(1
3, 1
3,1), and q|EAGF (2
3, 2
3,1) = q|EAGH(2
3, 2
3,1);
cf. Figure 8 and (3.27). To matchq at these two points, we let
(3.27)
v3,1 = 2
3q|EAGF (1
3, 1
3,1)
⎛
⎝
0
φ3,1
φ3,1 −
√
2φ3,3
⎞
⎠,
v3,2 = 1
3q|EAGF (2
3, 2
3,1)
⎛
⎝
0
−4φ3,1√
2φ3,3 −4φ3,1
⎞
⎠.
In order to match the other 8q values on the other two edges, we need to “borrow”
one degree of freedom ofvh from the next interface. Similar to (3.27), we construct
one more basis function on the next two tetrahedra:
(3.28) φ3,4 =
{
27
4 x(1−z)(y −x)(1−y)o n AGHE,
27
4 x(1−y)(z −x)(1−y)o n AGHD.
We only use one additional freedom, in addition to the 9-dimensional space (3.26)
⎛
⎝
0
φ3,4
0
⎞
⎠
whose divergence is zero on all edges except on theEAGF side of of edgeAG.W e
construct v3,i so that divv3,i match the rest of the eight degrees of freedom ofq at
the other two edges: edgeAG:
(3.29)
v3,3 = 2
3q|EAGF (1
3, 1
3, 1
3)
⎛
⎝
0
φ3,1 −φ3,4
(1/
√
2)φ3,2
⎞
⎠,
v3,4 = 2
3q|EAGH(1
3, 1
3, 1
3)
⎛
⎝
φ3,1
φ3,4
0
⎞
⎠,
v3,5 = −1
3q|EAGF (2
3, 2
3, 2
3)
⎛
⎝
0
4φ3,1 −φ3,4
(1/
√
2)φ3,2
⎞
⎠,
v3,6 = −1
3q|EAGH(2
3, 2
3, 2
3)
⎛
⎝
4φ3,1
φ3,4
0
⎞
⎠,

DIVERGENCE-FREE FINITE ELEMENTS ON TETRAHEDRAL GRIDS 683
and edgeAE:
(3.30)
v3,7 =
√
2
3 q|EAGF (0,0, 1
3)
⎛
⎝
φ3,2 −φ3,3√
2φ3,4
0
⎞
⎠,
v3,8 =
√
2
3 q|EAGH(0,0, 1
3)
⎛
⎝
0
φ3,2 −φ3,3 −
√
2φ3,4
φ3,2
⎞
⎠,
v3,9 = − 1
3
√
2
q|EAGF (0,0, 2
2)
⎛
⎝
φ3,2 −4φ3,3√
2φ3,4
0
⎞
⎠,
v3,10 = −3
√
2
q |EAGH(0,0, 2
3)
⎛
⎝
0
φ3,2 −4φ3,3 −
√
2φ3,4
φ3,2
⎞
⎠.
Hence, lettingv3 = ∑10
i=1 v3,i,w eh a v ed i vv3 zero at all vertices, and on all other
edges except the three edges on two sides of triangleAEG where the divergence
matches q, assumingEG is a boundary edge. We remark that we have to “borrow”
a degree of freedom from next internal triangle, no matter how high the polynomial
degree k is. For example, when k =5 ,w ed oh a v e3× 6 = 18 nodal degrees
of freedom forvh internal to three edges of triangleAEG, similar to (3.26) (see
Figure 9), whileq on the two sides ofAEG has 18 nodal values (recall that forP4,
we have dimensions 9 and 10 for them.) But{divvh} is still short of one dimension.
Now, for allk> 4, we have (k −4) mid-edge degrees of freedom on each edge,
shown in Figure 9. We ﬁrst construct a v3,m to match q values at the (k − 4)
mid-edge points. For example, fork =5 ,o ne d g eEA of Figure 8, we deﬁne
(3.31) φ3,m =
{
16y(z −x)2(1−z)2 on EAGH,
16x(z −y)2(1−z)2 on EAGF.
Then the divergence of
v3,m = q|EAGH(0,0, 1
2)
⎛
⎝
0
φ3,m
0
⎞
⎠
is zero on all edges except on the sideEAGH of edgeEA.A s
(q −divv3,m)EAGH(0,0, 1
2)=0 ,
we construct avh as (3.27)–(3.30) so that divvh|EAGH matches (q−divv3,m)EAGH
at two outside Lagrange nodes, (0,0, 1
4)a n d( 0,0, 3
4). For k> 5, we have exact
internal degrees of freedom for deﬁning (3.31) to matchq at internal edge nodes.
Hence, a construction can be done for edgesEA and AG for allk ≥ 4.
Next, as in the last lemma, we have to preserve the mean divergence-zero ele-
mentwise by correcting divv3 on three tetrahedra with twoP6 bubble functions,
supported on two neighboring tetrahedra each, as (3.9). For example,
b3,0 =
{
c0nAEGL2
AEHL2
EGHL2
GAH on EAGH,
nAEGL2
AEF L2
EGF L2
GAF on EAGF.

684 SHANGYOU ZHANG
E
F
B
A
G
H
D
C
E
F
B
A
G
H
D
C
Y
Z
W
X
Figure 10. The interior and boundary vertices of Ωh.
Here we need the divergence ofP6 bubble functions be zero at all vertices as well
as on all edges.
Repeating this construction for each of the 6 triangles around the diagonal edge
AG,w em a t c hq at all edges, assuming Ωh has only 6 tetrahedra. For a general
(small) cube ABCDEFGH in a reﬁned Ωh, the cube has 6 two-tetrahedra edges
like EA where two tetrahedra form a 90-degree face angle; cf. Figure 4. The
above construction would matchq exactly at 6 such two-tetrahedra edges. But the
(small) cube ABCDEFGH has also 6 one-tetrahedron edges likeEH and 6 ﬂat
two-tetrahedra edges likeEG (where two tetrahedra form a 180-degree face angle.)
For these 6 one-tetrahedron and 6 ﬂat two-tetrahedra edges, the above construction
may not match divv3,i with q there. We need to construct furtherv3,i for these
two cases.
If EH is a boundary edge, but inside a face square, such asQP in Figure 4, then
we use basis functions like (3.25), internal to triangleQPN to matchq|QP at the
bottom, to get av3,b.T h e n(q − divv3,b) would change theq values at the edge
QP on the two tetrahedra inside the top cube,QPNJ and QPJH . So we need to
repeat the work in (3.25)–(3.30) on the triangleQPJ.N e x t ,i fEH is an internal
edge, such asMN in Figure 4, theq values at edgeMN are matched separately
on the two cubes in front, and two cubes behind.
Finally, we consider the case of square-diagonal edge when it is not on the bound-
ary, for example,EG in Figure 4. There,AH and GD are such singular edges in
the other two directions. For simplicity of notation, we consider the case ofGD
depicted in Figure 10, where we assumeD is the origin, and the two cubes sharing
D are unit ones. HereGD is the intersection of two planes,ZGAD and HGCD.
We ﬁrst see why it is called a singular edge. At any pointx0 internal to the edge
GD, for anyvh ∈ Vh,k,w ew r i t e
(3.32) vh = u1
⎛
⎝
0
1/
√
2
1/
√
2
⎞
⎠+u2
⎛
⎝
1
0
0
⎞
⎠+u3
⎛
⎝
0
−1/
√
2
1/
√
2
⎞
⎠ =: u1 +u2 +u3,
where u1 is a globalPk polynomial on four tetrahedra sharing edgeGD, whileu2
and u3 are continuous piecewise-Pk on the four tetrahedra. By the continuity of

DIVERGENCE-FREE FINITE ELEMENTS ON TETRAHEDRAL GRIDS 685
vh,w eh a v e
divu1|AGDH(x0)=d i vu1|AGDC(x0)=d i vu1|ZGDH(x0)=d i vu1|ZGDC(x0),
divu2|AGDH(x0)=d i vu2|ZGDH(x0), divu2|AGDC(x0)=d i vu2|ZGDC(x0),
divu3|AGDH(x0)=d i vu3|AGDC(x0), divu3|ZGDH(x0)=d i vu3|ZGDC(x0).
Therefore, the dimension of the linear vector space of
(3.33) {divvh|AGDH(x0),divvh|AGDC(x0),divvh|ZGDH(x0),divvh|ZGDC(x0)}
is 3, not 4. Thus, for any point on edgeGD, we have a checkerboard mode, which
is limited by the constraint (cf. [11])
(3.34)
∑
i
(−1)i divvh|Ti(x0)=0 ,
where Ti stands for one of four tetrahedra around the edge GD. We need to
construct local basis functions for each of three linearly independent vectors in
(3.33). (3.32) provides a construction method. Let us consider ﬁrst theP4 case.
Let
φ3,s1 =
{
27
4 (x−y)(z −x)(1−z)(1−z), on ZGDH,
27
4 x(z −x)(1+ y −z)(1−z), on AGDH,(3.35)
φ3,s2 =
{
27
2 (x−y)(z −x)(1−z)x, on ZGDH,
27
2 x(z −x)(1+ y −z)x, on AGDH.(3.36)
We next deﬁne
v3,s1 = q|ZGDH (1
3,0, 1
3)1
3
⎛
⎝
0
0
4φ3,s1 −φ3,s2
⎞
⎠,
v3,s2 = q|ZGDH (2
3,0, 2
3)2
3
⎛
⎝
0
0
φ3,s2 −φ3,s1
⎞
⎠.
Then, div(v3,s1 + v3,s2)m a t c h e sq at the two Lagrange points on the edgeGD,
in ZGDH. Note that div(v3,s1 + v3,s2) = 0 at the two Lagrange points, on the
other side of planeAGDZ. Similarly, we can deﬁnev3,s3 and v3,s4 so that their
divergence matchesq at the two Lagrange points insideZGDH, while not altering
the match done on the other side of planeAGDZ. Hence
(3.37) q3 := q −div(v3,s1 +v3,s2 +v3,s3 +v3,s4)
vanishes on the edgeGD on y> 0 side. As we did for the non-singular edge case
AG, the construction (3.35)–(3.37) can be extended to anyPk, k ≥ 4. Again, we
correct q3 on each element to keep (3.21) byP6 bubble functions. By (3.34) and
(3.37), edge GD behaves as a boundary edge forq3. Hence the edge values ofq3
can be matched now by the divergence ofv3,i, deﬁned in (3.27), (3.29) and (3.30).
Summing over all suchv3,i over all edges of Ωh, after adding bubbles to preserve
(3.21), denoted byv3, it satisﬁes (3.20)–(3.22). □
After we match the element integrals, the vertex values and the edge values of
q ∈ Ph, we will next matchq on each face of element. This is the simplest task
among the others. The reason for this is that we can show the next lemma on any

686 SHANGYOU ZHANG
ˆK :
-
6
 
 
  	
           











@
@
@
@
@
@
@
@
@
ˆx
ˆy
ˆz
F(ˆx)
→
X X X X X X X X X X X X
       








Q Q Q Q Q Q Q Q Q
S
S
S
S
S
S
S
S
S
S
S
S
A
B
C
D
Figure 11. The aﬃne mapping fromˆK to K.
tetrahedral grid, unlike the other lemmas which are shown for the uniform grid
only.
Lemma 3.4. For anyq ∈ Ph deﬁned in(2.4) such that
∫
K q =0 ∀K ∈ Ωh, k ≥ 4
and q vanishes at all edges of gridΩh, there is a functionv4 ∈ Vh,k such that
divv4|TK
i
= q|TK
i
∀K ∈ Ωh,(3.38)
∫
K
divv4 =0 ∀K ∈ Ωh,(3.39)
‖v4‖H1(Ω)3 ≤ C‖q‖L2(Ω).(3.40)
Here TK
i , 1 ≤ i ≤ 4, are the four face triangles of tetrahedronK.
Proof. We note that fork =1 ,2,3, the lemma holds withv4 = 0 as q ≡ 0. To
understand the analysis better, we ﬁrst discuss the casek = 4, which is also covered
in the proof for generalk ≥ 4b e l o w .F o rk =4 ,l e tv4 = c0φK,w h e r eφK is the
bubble function of P4 on K.T h e n d i vv4 is a P3 function with zero integral on
K and zero trace on the 6 edges ofK. This is exactly howq is restricted in the
lemma. The three choices in c0 for v4 will provide a unique match to the three
degrees of freedom in deﬁningq on K. We next formalize this argument rigorously
for allk ≥ 4.
For a givenq speciﬁed in the lemma, we are going to constructv4 in 4 steps:
(3.41) v4 =( v4,1 +v4,2 +v4,3 +v4,4)φK ∈ C0(K)∩P3
k,
where v4,i are vector Pk−4 polynomials to be speciﬁed andφK is the P4 bubble
function on K.L e tK = ABCD with 4 face triangles numbered asT1 = ABC,
T2 = ABD, T3 = ACD and T4 = BCD.L e tF(ˆx)= Bˆx+x0 be an aﬃne mapping
from the referencetetrahedron ˆK = {0 ≤ ˆz ≤ 1−ˆx−ˆy, 0 ≤ y ≤ 1−ˆx, 0 ≤ ˆx ≤ 1}to
K so that the face triangles ofˆK on the plane ˆx =0 ,ˆy =0 ,ˆz =0a n dˆx+ˆy+ˆz =1
are mapped toT1, T2, T3 and T4, respectively. This is shown in Figure 11. Mapping
the equation div(v4,4)= q back to the reference element, we have
(3.42) ˆvT
4,4B−T∇φˆK +φˆKtrace
(
B−T∇ˆv4,4
)
=ˆq(ˆx),

DIVERGENCE-FREE FINITE ELEMENTS ON TETRAHEDRAL GRIDS 687
where
∇φˆK =
⎛
⎝
ˆyˆz(1−2ˆx− ˆy − ˆz)
ˆxˆz(1− ˆx−2ˆy − ˆz)
ˆxˆy(1− ˆx− ˆy −2ˆz)
⎞
⎠.
We choose, ifk> 4,
(3.43) ˆv4,4 = B
⎛
⎝
ˆx
ˆy
ˆz
⎞
⎠u4, for someu4 ∈ Pk−5.
Let ˆxi = ⟨ˆxi,ˆyi,ˆzi⟩be interior Lagrange points forPk−1 on the face triangleˆT4 =
F−1(T4) on the plane ˆx+ˆy+ˆz = 1 of the reference element. We derive the following
from (3.42):
(3.44) u4(ˆxi)= − ˆq(ˆxi)
ˆxiˆyiˆzi
.
Since ˆq = 0 on the three edges of the triangle and mapping back toK, we conclude
that
divv4,4|T4 = q|T4,
divv4,4|Ti =0 ,i ⁄=1 ,
|v4,4|H1 ≤ C‖q‖L2.
Now, fork = 4 in (3.43), we have to matchq on all four faces with only onev4,4
by letting
ˆv4,4 = Bc0, where c0 =
[
∇φˆK(ˆxi)
]−T
⎛
⎝
q(ˆx1)
q(ˆx2)
q(ˆx3)
⎞
⎠,
where ˆxi are the barycentric centers of any three-face triangle.
We repeat the construction ofv4,4 t h r e em o r et i m e st og e tv4,1, v4,2, v4,3 in
(3.41), whose divergence matchq on the other three triangles. Similar to (3.43) we
let
ˆv4,1 = B
⎛
⎝
1− ˆx− ˆy − ˆz
0
0
⎞
⎠u1, for someu1 ∈ Pk−5,
ˆv4,2 = B
⎛
⎝
0
1− ˆx− ˆy − ˆz
0
⎞
⎠u2, for someu2 ∈ Pk−5,
ˆv4,3 = B
⎛
⎝
0
0
1− ˆx− ˆy − ˆz
⎞
⎠u3, for someu3 ∈ Pk−5,

688 SHANGYOU ZHANG
where (cf. (3.44))ui are determined by the nodal values:
u1(ˆxi)= ˆq(ˆxi)
ˆyiˆzi(1− ˆyi − ˆzi)2 ∀ˆxi ∈ ˆT0
1 ,
u2(ˆxi)= ˆq(ˆxi)
ˆxiˆzi(1− ˆxi − ˆzi)2 ∀ˆxi ∈ ˆT0
2 ,
u3(ˆxi)= ˆq(ˆxi)
ˆxiˆyi(1− ˆxi − ˆyi)2 ∀ˆxi ∈ ˆT0
3 .
By the inverse reference mapping, we getv4,i. Lettingv4 = ∑4
i=1 v4,i. The lemma
is proven. □
Lemma 3.5. Let q ∈ Ph deﬁned in(2.4) with k ≥ 4 such that
∫
K q =0 ∀K ∈ Ωh
and q vanishes on all triangular faces of gridΩh. There is a functionv5 ∈ Vh,k
such that
divv5(x,y)= q(x,y) ∀(x,y) ∈ Ω,(3.45)
v5(x,y)=0 ∀(x,y) ∈ ∂K and ∀K ∈ Ωh,(3.46)
‖v5‖H1(Ω)3 ≤ C‖q‖L2(Ω).(3.47)
Proof. Each K ∈ Ωh is a scaling of one of 6 unit tetrahedra shown in Figure 4. The
properties listed in (3.45)–(3.47) are independent of scaling. Because we can work
out the other 5 cases similarly, we show one case in whichK is the unit tetrahedron
AEGH shown in Figure 4:
K = AEGH = {(x,y,z ) | 0 ≤ x ≤ y, 0 ≤ y ≤ z, 0 ≤ z ≤ 1}.
When restricted onK, we let the pressure space satisfying the lemma be
(3.48) PK =
{
q ∈ Pk−1(K)
⏐⏐⏐
∫
K
q =0 ,q |∂K =0
}
.
For anyq0 ∈ PK,w eh a v e
(3.49) q0 = φKq1, for someq1 ∈ Pk−5,
where φK is the degree-4 polynomial bubble function:
φK = λ1λ2λ3λ4, where(3.50)
λ1 = x, λ 2 = y −x, λ 3 = z −y,(3.51)
and λ4 =1 −λ1 −λ2 −λ3 =1 −z.
Further, q =0f o ra n yk ≤ 5, due to
∫
K q = 0, i.e., the lemma holds trivially for
4 ≤ k ≤ 5. We next introduce a subspace ofH1
0(K)3 ∩P3
k whose image under the
divergence operator is inside the spacePK deﬁned in (3.48):
(3.52) VK =
{
v ∈ Pk(K)3
⏐⏐⏐v|∂K =0 , divv ∈ PK
}
.
We show next that the divergence operator is also an onto mapping fromVK to
PK. Because of divergence-free polynomials, we would reduce the spaceVK to a
much smaller one which is mapped to the spacePK one-to-one by the divergence
operator. To make divv ∈ PK, we can limit
v = φK
⎛
⎝
λ1λ2v1(λ1,λ2,λ3)
λ2λ3v2(λ1,λ2,λ3)
λ3λ4v3(λ1,λ2,λ3)
⎞
⎠ where vi ∈ Pk−6(λ1,λ2,λ3),

DIVERGENCE-FREE FINITE ELEMENTS ON TETRAHEDRAL GRIDS 689
so that divv|∂K = 0. This can be seen by the calculation (3.54) below. After
eliminating divergence-free functions, we let
(3.53) V0 =
{
v ∈ Pk(K)3
⏐⏐⏐v = φK
⎛
⎝
λ1λ2v1(λ1,λ2,λ3)
λ2λ3v2(λ1,λ3)
λ3λ4v3(λ1)
⎞
⎠,v i ∈ Pk−6(K)4−i
}
.
Indeed, we have divv0 ∈ PK for v0 ∈ V0,a s
divv0 = φK[2(λ2 −λ1)v1 +λ1λ2v1x +2(1 −λ1 −2λ2 −λ3)v2(3.54)
+2v3(λ1 +λ2 +2λ3 −1)].
Now, for eachq0 ∈ PK, we will ﬁnd av0 ∈ V0 such that divv0 = q0.T h i si sd o n e
by mathematical induction. We ﬁrst construct av0 such that the highest order
terms of divv0 match those of a givenq0 ∈ PK. For anyq0 ∈ PK, we separate the
degree k −5t e r m so fq1 in (3.48) from the rest as follows:
(3.55) q1 =
k−5∑
i=0
k−5−i∑
j=0
∑
l=k−5−i−j
qijlλi
1λj
2λl
3 +q2(λ1,λ2,λ3),
where q2 is a degree (k−6) polynomial. When we compare the degreek−5t e r m s
of divv0 in (3.54) andq1 in (3.55), we need to check only the degree (k−6) terms
in v1, v2 and v3.W el e tav0 in (3.53) be
v1(λ1,λ2,λ3)=
k−6∑
i=0
k−6−i∑
j=0
∑
l=k−6−i−j
v1,ijl λi
1λj
2λl
3,
v2(λ1,λ3)=
k−6∑
i=0
∑
l=k−6−i
v2,i0l λi
1λl
3,(3.56)
v3(λ1)=
∑
i=k−6
v3,i00 λi
1.
We note that there are
k−5∑
i=0
k−5−i∑
j=0
1=
k−5∑
i=0
(i+1)= (k −3)(k−4)
2
coeﬃcients ofqijl in (3.55), which deﬁnes the (k −3)(k−4)/2 linear equations for
the unknown coeﬃcients ofvi in (3.56):
(
k−6∑
i=0
k−6−i∑
j=0
1)+(
k−6∑
i=0
1)+1 = (k −4)(k −5)
2 +k −5+1= (k −3)(k −4)
2 .
Wecanlist the( k−3)(k−4)/2linearequationsin thefollowing ordertogetan upper
triangular system except the last three equations involvingv1,(k−6)00, v2,(k−6)00 and
v3,(k−6)00:
For i =0 ,
−2v2,00(l−1) = qijl,j =0 ,(3.57)
2v1,00l = qijl +4v2,00l,j =1 ,(3.58)
2v1,0(j−1)l = qijl,j =2:( k −5).(3.59)

690 SHANGYOU ZHANG
For i =1:( k −7),
−2v2,i0(l−1) = qijl +2v1,(i−1)0l −2v2,(i−1)0l,j =0 ,(3.60)
(2+ i)v1,i0l = qijl +2v1,(i−1)1l +4v2,i0l,j =1 ,(3.61)
(2+ i)v1,i(j−1)l = qijl +2v1,(i−1)jl,j =2:( k −5−i).(3.62)
For i = k −6,
−2v2,i00 +4v3,i00 = qijl +2v1,(i−1)0l −2v2,(i−1)0l,j =0 ,(3.63)
(2+ i)v1,i00 +2v3,i00 = qijl +2v1,(i−1)1l +4v2,i0l,j =1 .(3.64)
Finally, fori = k −5,
−2v1,(i+1)00 −2v2,(i+1)00 +v2,(i+1)00 = qijl +2v1,(i−1)jl.(3.65)
Here in equations (3.57)–(3.65), the indexl = k − 5 − i − j. Also in all of these
equations, allvm,ijl on the right-hand side are resolved by earlier equations. For the
last three equations (3.63)–(3.65), the determinant of the 3× 3 coeﬃcient matrix
is −4(k −6). So, for ak ≥ 7, we ﬁnd a uniquev0 so that divv0 and q0 match the
highest order terms. We move the divv0 to the right-hand side of (3.48), combined
into q2 there. We then repeat the above construction for one lower degreeq0,u n t i l
k =6 .W h e nk = 6, the systems of equations (3.57)–(3.65) become to
−2v2,000 = q001,
2v1,000 = q010 +4v2,000,
2v3,000 = q100 +2v1,000 +2v2,000.
This system is an upper triangular one, and has a unique solution.
Therefore, we constructed a locally supportedv5 on one tetrahedron AEGH.
Similar construction can be done on the other ﬁve types of tetrahedra. The lemma
is proven. □
Corollary 3.1. Let k ≥ 6. The mixed ﬁnite element pair(Vh,k,Ph) are deﬁned
in (2.3) and (2.4).L e tn be the number of cubes in each coordinate direction; cf.
Figure 2. The dimensions ofVh,k, Ph and subspaceZh (deﬁned in(2.8)) are
dimVh,k =( nk −1)3,(3.66)
dimPh = 1
6n3(k +2)(k +1)k −3kn(n2 +n+2)+5 ,(3.67)
dimZh =d i mVh,k −dimPh.(3.68)

DIVERGENCE-FREE FINITE ELEMENTS ON TETRAHEDRAL GRIDS 691
Proof. We summarize all constraints for divwh = q in the last few lemmas:
Corner vertices – Type (b): 2 ×6,
Corner vertices – Type (a): 3 ×2,
Mid-edge vertices – Type (c): 4 ×6(n−1),
Mid-edge vertices – Type (d): 3 ×6(n−1),
Mid-face vertices – Type (e): 4 ×6(n−1)2,
Internal vertices – Type (f): 6 ×(n−1)3,
One-tetrahedron boundary edges: ( k −2)×6n,
Diagonal boundary edges: ( k −2)×6n2,
Internal singular edges: ( k −2)×3(n−1)n2,
Global integral constraint: 1 .
Deducting the number of constraints from the dimension of discontinuousPk−1
polynomials on Ωh,w ep r o v et h ec o r o l l a r y . □
Numerically, we have veriﬁed Corollary 3.1:
dimPh =
⎧
⎪⎪⎪⎨
⎪⎪⎪⎩
269 if k =6 ,n =1 ,
2405 if k =6 ,n =2 ,
425 if k =7 ,n =1 ,
3701 if k =7 ,n =2 .
We proved (3.67) fork ≥ 6 only, but it seems to hold fork = 5 as well:
dimPh = 155 and 1445, if n =1a n d2, respectively.
However, (3.67) no longer holds fork ≤ 4:
dimPh = 75 and 772, if n =1a n d2, respectively,
while (3.67) gives 76 and 789, respectively.
Theorem 3.1. Let k ≥ 6. The mixed ﬁnite element pair (Vh,k,Ph) deﬁned in
(2.3) and (2.4) is stable on the uniform grids, i.e., the following inf-sup condition
holds:
(3.69) inf
q⁄=0,q∈Ph
sup
vh∈Vh,k
b(vh,q)
‖vh‖H1(Ω)3‖q‖L2(Ω)
≥ C.
Proof. For anyq ∈ Ph,w ec o n s t r u c tavh ∈ Vh,k to satisfy (3.69). By (3.1), there
is av1 ∈ Vh such that
∫
K
(q −divv1)=0 ∀K ∈ Ωh,
‖v1‖H1 ≤ C1‖q‖L2.
By (3.2), there isv2 ∈ Vh such that
[divv2 −q +div v1]|K (aK
i )=0 ∀K ∈ Ωh and for all vertices ofK,
‖v2‖H1 ≤ C2‖q −divv1‖L2.

692 SHANGYOU ZHANG
By (3.20), there isv3 ∈ Vh such that
[divv3 −q +div(v1 +v2)]|K (EK
i )=0 ∀K ∈ Ωh and for all edges ofK,
‖v3‖H1 ≤ C3‖q −div(v1 −v2)‖L2.
By (3.38), there isv4 ∈ Vh such that
[divv4 −q +div(v1 +v2 +v3)]|K(FK
i )=0 ∀K ∈ Ωh
and for all face triangles ofK,
‖v4‖H1 ≤ C4‖q −div(v1 +v2 +v3)‖L2.
By (3.45), there is av5 ∈ Vh such that
divv5 = q −div(v1 +v2 +v3 +v4),
‖v5‖H1 ≤ C5‖q −div(v1 +v2 +v3 +v4)‖L2.
Let v = −v1 −v2 −v3 −v4 −v5. It follows that
‖v‖H1 ≤‖v1‖H1 +‖v2‖H1 +‖v3‖H1 +‖v4‖H1 +‖v5‖H1
≤ C1‖q‖L2 +C2‖q −divv1‖L2 +C3‖q −div(v1 +v2)‖L2 +···
≤ C1‖q‖L2 +C2(‖q‖L2 +‖divv1‖L2)+ ···
≤ C1‖q‖L2 +C2(‖q‖L2 +C1‖q‖L2)+ ···
≤ C∗‖q‖L2
and that
b(vh,q)=( −divvh,q)= ‖q‖2
L2(Ω) ≥ C−1
∗ ‖v‖H1(Ω)3‖q‖L2(Ω).
(3.69) is proved withC = C−1
∗ . □
Theorem 3.2. Let k ≥ 6. The discrete solution(uh,ph) of (2.5) approximate that
of (2.2) in the optimal order:
‖u−uh‖H1(Ω)3 +‖p−ph‖L2(Ω)(3.70)
≤ Chmin{k,r}(‖u‖Hr+1(Ω)3 +‖p‖Hr(Ω)),r ≥ 1.
Proof. By the inf-sup condition (3.69) and the standard mixed ﬁnite element theory
[10], it follows that
‖u−uh‖H1(Ω)3 +‖p−ph‖L2(Ω)
≤ C(i n f
vh∈Vh,k
‖u−vh‖H1(Ω)3 +i n f
qh∈Ph
‖p−qh‖L2(Ω))
≤ C(i n f
vh∈Vh,k
‖u−vh‖H1(Ω)3 +i n f
qh∈˜Ph
‖p−qh‖L2(Ω))
where ˜Ph is the space of continuousPk−1 polynomials with mean value zero:
˜Ph = Ph ∩C(Ω) =
{
qh ∈ C(Ω) |
∫
Ω
qh =0 ,q h|K ∈ Pk−1 ∀K ∈ Ωh
}
.
The theorem is proven as both spacesVh,k and ˜Ph provide the optimal order of
approximation. □

DIVERGENCE-FREE FINITE ELEMENTS ON TETRAHEDRAL GRIDS 693
4. Numerical tests
In this section, we report some numerical tests on thePk-Pk−1 elements for the
stationary Stokes equations (2.1) on the unit cube, Ω = (0,1)3.T h e g r i d s a r e
obtained by the standard multigrid reﬁnement; cf. [16]. The ﬁrst three grids are
depicted in Figure 2.
Figure 12. The exact solution, the ﬁrst component ofu and p
in (4.2), restricted onz =0 .33.
We choose the right-hand side functionf for (2.1) as
(4.1)
f = −Δcurl
⎛
⎝
0
g
g
⎞
⎠+ 1
9∇gxy
=
⎛
⎝
−gxxy −gyyy −gyzz +gxxz +gyyz +gzzz +gxxy/9
−gxxx −gxyy −gxzz +gxyy/9
gxxx +gxyy +gxzz +gxyz/9
⎞
⎠,
where
g =2 12(x−x2)2(y −y2)2(z −z2)2.
The exact solution for the Stokes equations (2.1) is
(4.2) u = curl
⎛
⎝
0
g
g
⎞
⎠,p = 1
9gxy.
As we are unable to plot a 3D function in 4D, we show the restriction of the
functions u (the ﬁrst component) andp, on the planez =0 .33 in Figure 12. We
note that the grids obtained by the intersection of tetrahedra in Ωh and the plane
consist of both rectangles and triangles, shown at the bottom in Figure 12 and in
Figure 13.
In Table 1 we list errors for thePk-Pk−1 element fork = 6, on three level of grids
Ωh. The iterated penalty method is used to solve the discrete linear equations. The
order of convergence ﬁts the estimate (3.70) well. We show some errors in Figure 14.
Table 1. The errors for thePk-Pk−1 ( k = 6) element on Figure 2 grids.
|u−uh|H1 hn ‖p−ph‖L2 hn
1 6.73310 29.66007
2 0.23981 4.81 1.13377 4.70
3 0.00421 5.83 0.02196 5.69

694 SHANGYOU ZHANG
Figure 13. T h ec u to nt h et h i r dl e v e lg r i dΩh by planez =0 .33.
Figure 14. The errors for the ﬁrst component of u and p re-
stricted on planez =0 .33.
Acknowledgments
This work was initially supported by the National Science Foundation Award
9625907 and ﬁnished in June 2008 during the author’s visit to Dr. Xuejun Xu,
sponsored by the State Key Laboratory of Scientiﬁc and Engineering Computing,
Beijing, China.
The author thanks an anonymous referee who pointed out numerous mistakes in
an early version of this manuscript.
References
[1] D. N. Arnold and J. Qin,Quadratic velocity/linear pressure Stokes elements, in Advances in
Computer Methods for Partial Diﬀerential Equations VII, R. Vichnevetsky and R.S. Steple-
men, eds., 1992.
[2] D. Boﬃ,Three-dimensional ﬁnite element methods for the Stokes problem,SIAM J. Numer.
Anal. 34 (1997), 664–670. MR1442933 (98a:65160)
[3] S. C. Brenner and L. R. Scott, The Mathematical Theory of Finite Element Methods,
Springer-Verlag, New York, 1994. MR1278258 (95f:65001)
[4] F. Brezzi and M. Fortin, Mixed and hybrid ﬁnite element methods, Springer, 1991.
MR1115205 (92d:65187)
[5] P.G. Ciarlet, The Finite ElementMethod forElliptic Problems, North-Holland, Amsterdam,
1978. MR0520174 (58:25001)
[6] M. Fortin and R. Glowinski, Augmented Lagrangian Methods: Applications to the Numer-
ical Solution of Boundary-value Problems, North Holland, Amsterdam, 1983. MR724072
(85a:49004)

DIVERGENCE-FREE FINITE ELEMENTS ON TETRAHEDRAL GRIDS 695
[7] J. Pitk¨aranta and R. Stenberg,Error bounds for the approximation of the Stokes problem
using bilinear/constant elements on irregular quadrilateral meshes,i nT h eM a t h e m a t i c so f
ﬁnite elements and applications V, J. Whiteman, ed., Academic Press, London, 1985, 325–
334. MR811045
[8] J. Qin, On the convergence of some low order mixed ﬁnite elements for incompressible ﬂuids,
Thesis, Pennsylvania State University, 1994.
[9] J. Qin and S. Zhang,Stability and approximability of theP1-P0 element for Stokes equations,
Int. J. Numer. Meth. Fluids54 (2007), no. 5, 497–515. MR2322456 (2008b:65153)
[10] P. A. Raviart and V. Girault, Finite element methods for Navier-Stokes equations, Springer,
1986. MR851383 (88b:65129)
[11] L. R. Scott and M. Vogelius,Norm estimates for a maximal right inverse of the divergence
operator in spaces of piecewise polynomials, RAIRO, Modelisation Math. Anal. Numer. 19
(1985), 111–143. MR813691 (87i:65190)
[12] L. R. ScottandM. Vogelius,Conforming ﬁnite element methods for incompressible and nearly
incompressible continua, in Lectures in Applied Mathematics 22, 1985, 221–244. MR818790
(87h:65202)
[13] L. R. Scott and S. Zhang, Finite element interpolation of nonsmooth functions satisfying
boundary conditions, Math. Comp. 54 (1990), 483–493. MR1011446 (90j:65021)
[14] L. R. Scott and S. Zhang, Multilevel Iterated Penalty Method for Mixed Elements, the
Proceedings for the Ninth International Conference on Domain Decomposition Methods, 133-
139, Bergen, 1998.
[15] M. Vogelius, A right-inverse for the divergence operator in spaces of piecewise polynomials
application to the p version of the ﬁnite element method, Numer. Math. 41 (1983), 19–37.
MR696548 (85f:65113a)
[16] S. Zhang,Successive subdivisions of tetrahedra andmultigrid methods on tetrahedral meshes,
Houston J. of Math. 21 (1995), 541–556. MR1352605 (96f:65183)
[17] S. Zhang,A new family of stable mixed ﬁnite elements for 3D Stokes equations,M a t h .C o m p .
74 (2005), 250, 543–554. MR2114637 (2005j:65151)
[18] S. Zhang, On the P1 Powell-Sabin divergence-free ﬁnite element for the Stokes equations,
J. Comp. Math., 26 (2008), 456-470. MR2421893 (2009j:76160)
Department of Mathematical Sciences, University of Delaware, Newark, Delaware
19716
E-mail address: szhang@udel.edu