Ctrl+K

Finite Elements in H(\operatorname{div})

\(\DeclareMathOperator{\opdiv}{div}\)

48. Finite Elements in \(H(\operatorname{div})\)#

Finite element sub-spaces of \(H(\opdiv)\) are built from element-wise polynomials, which normal-components are continuous across element boundaries. This is obtained by proper degrees of freedom.

We introduce two families of spaces: The Raviart-Thomas (RT) elements, and the Brezzi-Douglas-Marini (BDM) elements.

The lowest order triangular element is the \(RT_0\) element. Here, the element space is

\[ V_T = \{ \vec a + b \vec x : \vec a \in {\mathbb R}^2, b \in {\mathbb R} \} \]

Its dimension is 3. Thus we need three degrees of freedom (what are linear functionals) to define a particular function from \(V_T\). These functionals are the integrals of the normal component over edges: $\( \tau \mapsto \int_E \tau_n \)$

Since \(\tau_n\) is constant on an edge, this value uniquely defines the normal component on the edge.

The next element is the \(BDM_1\) element: Its element space is

\[ V_T = [P^1]^2 \]

It has dimension 6. We need two functionals per edge to define the normal component.

The \(BDM_k\) elements are defined as

\[ V_T = [P^k]^2 \]

the degrees of freedom are

  • \(\int_E \tau_n q_i \quad \quad \) with \(q_i\) basis for \(P^k(E)\)

  • \(\int_T \opdiv \tau \, r_j \quad\) with \(r_j\) basis for \(P^{k-1}(T)/{\mathbb R}\)

  • \(\int_T \tau \operatorname{curl} s_l \quad \quad\) with \(s_j\) basis for \(P_0^{k+1}\)

The divergence of \(BDM_k\) elements give

\[ \opdiv V_T = P^{k-1} \]

Raviart-Thomas elements \(RT_k\) are similar

\[ V_T = \{ \vec a + b \vec x : \vec a \in [P^k]^2, b \in P^k \}. \]

There is

\[ [P^k]^2 \subset RT_k \subset [P^{k+1}]^2 \]

and

\[ \opdiv RT_k = P^k \]

The Raviart-Thomas space leading to the same range of the divergence is smaller than the BDM space.

48.1. Piola Transformation#

Typically, one defines a reference element, and maps it to the physical elements in the mesh. Basis functions are implemented for the reference element.

Mapping and its Jacobian:

\[ x = \Phi (\widehat x) \qquad F = \Phi^\prime \qquad J = \det F \]

The standard mapping for \(H^1\) elements is

\[ w(x) := \widehat w (\widehat x) \]

The gradient on the physical element is computed from the chain-rule:

\[ \nabla w (x) = F^{-T} \nabla \widehat w (\widehat x) \]

\(H(\opdiv)\) elements are mapped by the Piola transform, i.e.

\[ \tau(x) := J^{-1} F \, \widehat \tau (\widehat x) \]

From that, the divergence follows as algebraic transformation of the divergence on the reference element:

\[ \opdiv \tau (x) = J^{-1} \widehat \opdiv \widehat \tau (\widehat x) \]

Proof: technical calculation in strong from, short proof using the weak form of the divergence.

Furthermore, fluxes through edges are preserved:

\[ \int_E \tau_n = \int_{\widehat E} \widehat \tau_n \qquad \text{for} \quad E = \Phi (\widehat E) \]

Theorem: BDM/RT finite elements are interpolation equivalent when using the Piola transform, i.e.

\[ I_T ( J^{-1} F \widehat \sigma ) = J^{-1} F \, ( I_\widehat T \widehat \sigma ) \]
from ngsolve import *
from ngsolve.webgui import Draw
mesh = Mesh(unit_square.GenerateMesh(maxh=0.3))

# a RT0-space:
V = HDiv(mesh, order=0, RT=True)
print ("num edges =", mesh.nedge)
print ("ndof = ", V.ndof)
tau = GridFunction(V)
tau.vec[:] = 0 
tau.vec[35] = 1
Draw(tau, mesh, "tau", vectors = { "grid_size" : 30 });
Draw(div(tau), mesh);

num edges = 42
ndof =  42

High order RT and BDM spaces are built by an hierarchical basis. The contain the lowest order \(RT_0\) basis functions. Then, basis functions on edges, faces (and cells) are added.

See Dissertation Sabine Zaglmayr, page 81 for the triangle.

# a BDM2-space:
V = HDiv(mesh, order=2)
tau = GridFunction(V)

Except for the first one, all shape functions on edges are div-free:

print ("dofs on edge 22:",V.GetDofNrs(NodeId(EDGE,22)))
tau.vec[:] = 0
tau.vec[V.GetDofNrs(NodeId(EDGE,22))[1]] = 1

Draw(tau, mesh, vectors = { "grid_size" : 30 })
Draw(div(tau), mesh);

dofs on edge 22: (22, 86, 87)

print ("dofs on face 0:", V.GetDofNrs(NodeId(FACE,0)))

tau.vec[:] = 0
tau.vec[V.GetDofNrs(NodeId(FACE,0))[1]] = 1
Draw(tau, mesh, vectors = { "grid_size" : 30 })
Draw(div(tau), mesh);

dofs on face 0: (126, 127, 128)

With RT=True we get the Raviart-Thomas space, which is \([P^k]^d + P^k x\):

V = HDiv(mesh, order=2, RT=True)
tau = GridFunction(V)
func = y*x*CF( (x,y) )
tau.Set(func)
print ("err =", sqrt(Integrate((tau-func)**2, mesh)))

err = 4.723697328914656e-16
Contents