Non-conforming Finite Element Methods

25. Non-conforming Finite Element Methods#

In a conforming finite element method, one chooses a sub-space $V_{h} \subset V$ , and defines the finite element approximation as

Find u_{h} \in V_{h} : A (u_{h}, v_{h}) = f (v_{h}) \forall v_{h} \in V_{h}

For reasons of simpler implementation, or even of higher accuracy, the conforming framework is often violated. Examples are:

The finite element space $V_{h}$ is not a sub-space of $V = H^{m}$ . Examples are the non-conforming $P^{1}$ triangle, and the Morley element for approximation of $H^{2}$ .
The Dirichlet boundary conditions are interpolated in the boundary vertices.
The curved domain is approximated by straight sided elements
The bilinear-form and the linear-form are approximated by inexact numerical integration

The lemmas by Strang are the extension of Cea’s lemma to the non-conforming setting.

Ref: Gil Strang: Variational crimes in the finite element method

25.1. The First Lemma of Strang#

In the first step, let $V_{h} \subset V$ , but the bilinear-form and the linear-form are replaced by mesh-dependent forms

A_{h} (., .) : V_{h} \times V_{h} \to R

and

f_{h} (.) : V_{h} \to R .

We do not assume that $A_{h}$ and $f_{h}$ are defined on $V$ . We assume that the bilinear-forms $A_{h}$ are uniformly coercive, i.e., there exists an $α_{1}$ independent of the mesh-size such that

A_{h} (v_{h}, v_{h}) \geq α_{1} ‖ v_{h} ‖_{V}^{2} \forall v_{h} \in V_{h}

The finite element problem is defined as

Find u_{h} \in V_{h} : A_{h} (u_{h}, v_{h}) = f_{h} (v_{h}) \forall v_{h} \in V_{h}

First Lemma of Strang: Assume that

$A (., .)$ is continuous on $V$

$A_{h} (., .)$ is uniformly coercive

Then there holds

$\begin{array}{rcl} ‖ u - u_{h} ‖ & ⪯ & inf_{v_{h} \in V_{h}} {‖ u - v_{h} ‖ + sup_{w_{h} \in V_{h}} \frac{| A (v_{h}, w_{h}) - A_{h} (v_{h}, w_{h}) |}{‖ w_{h} ‖}} \\ + sup_{w_{h} \in V_{h}} \frac{f (w_{h}) - f_{h} (w_{h})}{‖ w_{h} ‖} \end{array}$

Proof: Choose an arbitrary $v_{h} \in V_{h}$ , and set $w_{h} := u_{h} - v_{h}$ . We use the uniform coercivity, and the definitions of $u$ and $u_{h}$ :

\begin{aligned} α_{1} ‖ u_{h} - v_{h} ‖_{V}^{2} & \leq A_{h} (u_{h} - v_{h}, u_{h} - v_{h}) = A_{h} (u_{h} - v_{h}, w_{h}) \\ = A (u - v_{h}, w_{h}) + [A (v_{h}, w_{h}) - A_{h} (v_{h}, w_{h})] + [A_{h} (u_{h}, w_{h}) - A (u, w_{h})] \\ = A (u - v_{h}, w_{h}) + [A (v_{h}, w_{h}) - A_{h} (v_{h}, w_{h})] + [f_{h} (w_{h}) - f (w_{h})] \end{aligned}

Divide by $‖ u_{h} - v_{h} ‖ = ‖ w_{h} ‖$ , and use the continuity of $A (., .)$ :

‖ u_{h} - v_{h} ‖ ⪯ ‖ u - v_{h} ‖ + \frac{| A (v_{h}, w_{h}) - A_{h} (v_{h}, w_{h}) |}{‖ w_{h} ‖} + \frac{| f (w_{h}) - f_{h} (w_{h}) |}{‖ w_{h} ‖}

Using the triangle inequality, the error $‖ u - u_{h} ‖$ is bounded by

\begin{aligned} ‖ u - u_{h} ‖ & \leq inf_{v_{h} \in V_{h}} ‖ u - v_{h} ‖ + ‖ v_{h} - u_{h} ‖ \\ ⪯ inf_{v_{h} \in V_{h}} {‖ u - v_{h} ‖ + \frac{| A (v_{h}, w_{h}) - A_{h} (v_{h}, w_{h}) |}{‖ w_{h} ‖} + \frac{| f (w_{h}) - f_{h} (w_{h}) |}{‖ w_{h} ‖}} \\ ⪯ inf_{v_{h} \in V_{h}} {‖ u - v_{h} ‖ + sup_{w_{h} \in V_{h}} \frac{| A (v_{h}, w_{h}) - A_{h} (v_{h}, w_{h}) |}{‖ w_{h} ‖} + sup_{w_{h} \in V_{h}} \frac{| f (w_{h}) - f_{h} (w_{h}) |}{‖ w_{h} ‖}} \end{aligned}

The lemma is proven. $◻$

Example: Lumping of the $L_{2}$ bilinear-form: Define the $H^{1}$ - bilinear-form

A (u, v) = \int_{Ω} \nabla u \cdot \nabla v + \int_{Ω} u v d x,

and perform Galerkin discretization with $P^{1}$ triangles. The second term leads to a non-diagonal matrix. The vertex integration rule

\int_{T} v d x \approx \frac{| T |}{3} \sum_{α = 1}^{3} v (x_{T, α})

is exact for $v \in P^{1}$ . We apply this integration rule for the term $\int u v d x$ :

A_{h} (u, v) = \int \nabla u \cdot \nabla v + \sum_{T \in T} \frac{| T |}{3} \sum_{α = 1}^{3} u (x_{T, α}) v (x_{T, α})

The bilinear form is now defined only for $u, v \in V_{h}$ , since point evaluation is not defined on $H^{1}$ for $d \geq 2$ . The integration is not exact, since $u v \in P^{2}$ on each triangle.

Inserting the nodal basis $φ_{i}$ , we obtain a diagonal matrix for the second term:

\begin{array}{r} φ_{i} (x_{T, α}) φ_{j} (x_{T, α}) = {\begin{array}{cl} 1 & for x_{i} = x_{j} = x_{T, α} \\ 0 & else \end{array} \end{array}

To apply the first lemma of Strang, we have to verify the uniform coercivity

\sum_{T} \frac{| T |}{3} \sum_{α = 1}^{3} | v_{h} (x_{T, α}) |^{2} \geq α_{1} \sum_{T} \int_{T} | v_{h} |^{2} d x \forall v_{h} \in V_{h},

which is done by transformation to the reference element. The consistency error can be estimated by (exercise!)

| \int_{T} u_{h} v_{h} d x - \frac{| T |}{3} \sum_{α = 1}^{3} u_{h} (x_{α}) v_{h} (x_{α}) | ⪯ h_{T}^{2} ‖ \nabla u_{h} ‖_{L_{2} (T)} ‖ \nabla v_{h} ‖_{L_{2} (T)}

Summation over the elements give

A (u_{h}, v_{h}) - A_{h} (u_{h}, v_{h}) ⪯ h^{2} ‖ u_{h} ‖_{H^{1} (Ω)} ‖ v_{h} ‖_{H^{1} (Ω)}

The first lemma of Strang proves that this modification of the bilinear-form preserves the order of the discretization error:

\begin{aligned} ‖ u - u_{h} ‖_{H^{1}} & ⪯ inf_{v_{h} \in V_{h}} {‖ u - v_{h} ‖_{H^{1}} + sup_{w_{h} \in V_{h}} \frac{| A (v_{h}, w_{h}) - A_{h} (v_{h}, w_{h}) |}{‖ w_{h} ‖_{H^{1}}}} \\ ⪯ ‖ u - I_{h} u ‖_{H^{1}} + sup_{w_{h} \in V_{h}} \frac{| A (I_{h} u, w_{h}) - A_{h} (I_{h} u, w_{h}) |}{‖ w_{h} ‖_{H^{1}}} \\ ⪯ h ‖ u ‖_{H^{2}} + sup_{w_{h} \in V_{h}} \frac{h^{2} ‖ I_{h} u ‖_{H^{1}} ‖ w_{h} ‖_{H^{1}}}{‖ w_{h} ‖_{H^{1}}} \\ ⪯ h ‖ u ‖_{H^{2}} \end{aligned}

A diagonal $L_{2}$ matrix has some advantages:

It avoids oscillations in boundary layers (exercises!)
In explicit time integration methods for parabolic or hyperbolic problems, one has to solve linear equations with the $L_{2}$ -matrix. This becomes cheap for diagonal matrices. See Mass lumping and local time-stepping

25.2. The Second Lemma of Strang#

In the following, we will also skip the requirement $V_{h} \subset V$ . Thus, the norm $‖ . ‖_{V}$ cannot be used on $V_{h}$ , and it will be replaced by mesh-dependent norms $‖ . ‖_{h}$ . These norms must be defined for $V + V_{h}$ . As well, the mesh-dependent forms $A_{h} (., .)$ and $f_{h} (.)$ are defined on $V + V_{h}$ . We assume

uniform coercivity:

$A_{h} (v_{h}, v_{h}) \geq α_{1} ‖ v_{h} ‖_{h}^{2} \forall v_{h} \in V_{h}$
continuity:

$A_{h} (u, v_{h}) \leq α_{2} ‖ u ‖_{h} ‖ v_{h} ‖_{h} \forall u \in V + V_{h}, \forall v_{h} \in V_{h}$

The error can now be measured only in the discrete norm $‖ u - u_{h} ‖_{V_{h}}$ .

Second Lemma of Strang:

Under the above assumptions there holds

$‖ u - u_{h} ‖_{h} ⪯ inf_{v_{h} \in V_{h}} ‖ u - v_{h} ‖_{h} + sup_{w_{h} \in V_{h}} \frac{| A_{h} (u, w_{h}) - f_{h} (w_{h}) |}{‖ w_{h} ‖_{h}}$

Remark: The first term in the estimate is the approximation error, the second one is called consistency error.

Proof: Let $v_{h} \in V_{h}$ . Again, set $w_{h} = u_{h} - v_{h}$ , and use the $V_{h}$ -coercivity:

\begin{aligned} α_{1} ‖ u_{h} - v_{h} ‖_{h}^{2} & \leq A_{h} (u_{h} - v_{h}, u_{h} - v_{h}) = A_{h} (u_{h} - v_{h}, w_{h}) \\ = A_{h} (u - v_{h}, w_{h}) + [f_{h} (w_{h}) - A_{h} (u, w_{h})] \end{aligned}

Again, divide by $‖ u_{h} - v_{h} ‖$ , and use continuity of $A_{h} (., .)$ :

‖ u_{h} - v_{h} ‖_{h} ⪯ ‖ u - v_{h} ‖_{h} + \frac{A_{h} (u, w_{h}) - f_{h} (w_{h})}{‖ w_{h} ‖_{h}}

The rest follows from the triangle inequality as in the first Lemma of Strang. $◻$

The non-conforming $P^{1}$ triangle

The non-conforming $P^{1}$ triangle is also called the Crouzeix-Raviart element.

The finite element space generated by the non-conforming $P^{1}$ element is

V_{h}^{n c} := {v \in L_{2} : v_{| T} \in P^{1} (T), and v is continuous in edge mid-points}

The functions in $V_{h}^{n c}$ are not continuous across edges, and thus, $V_{h}^{n c}$ is not a sub-space of $H^{1}$ . We have to extend the bilinear-form and the norm in the following way:

A_{h} (u, v) = \sum_{T \in T} \int_{T} \nabla u \nabla v d x \forall u, v \in V + V_{h}^{n c}

and

‖ v ‖_{h}^{2} := \sum_{T \in T} ‖ \nabla v ‖_{L_{2} (T)}^{2} \forall v \in V + V_{h}^{n c}

We consider the Dirichlet-problem with $u = 0$ on $Γ_{D}$ . Thus, $‖ \cdot ‖_{V_{h}}$ is a norm.

from ngsolve import *
from ngsolve.webgui import Draw

mesh = Mesh(unit_square.GenerateMesh(maxh=0.15))

fes = FESpace("nonconforming", mesh, order=1, dirichlet=".*")
u,v = fes.TnT()

a = BilinearForm(grad(u)*grad(v)*dx).Assemble()
f = LinearForm(5*v*dx).Assemble()

gfu = GridFunction(fes)
gfu.vec.data = a.mat.Inverse(freedofs=fes.FreeDofs())*f.vec
Draw (gfu, deformation=True);

We will apply the second lemma of Strang.

The continuous $P^{1}$ finite element space $V_{h}^{c}$ is a sub-space of $V_{h}^{n c}$ . Let $I_{h} : H^{2} \to V_{h}^{c}$ be the nodal interpolation operator.

To bound the approximation term in Strang’s second lemma, we use the inclusion $V_{h}^{c} \subset V_{h}^{n c}$ :

inf_{v_{h} \in V_{h}^{n c}} ‖ u - v_{h} ‖_{h} \leq ‖ u - I_{h} u ‖_{H^{1}} ⪯ h ‖ u ‖_{H^{2}}

To bound the consistency error, we perform integration by parts on every element:

\begin{array}{rcl} r (w_{h}) & = & A_{h} (u, w_{h}) - f (w_{h}) \\ = & \sum_{T} \int_{T} \nabla u \nabla w_{h} - \sum_{T} \int_{T} f w_{h} d x \\ = & \sum_{T} \int_{\partial T} \frac{\partial u}{\partial n} w_{h} d s - \sum_{T} \int_{T} (Δ u + f) w_{h} d s \\ = & \sum_{T} \int_{\partial T} \frac{\partial u}{\partial n} w_{h} d s \end{array}

Let $E$ be an edge of the triangle $T$ . Define the mean value ${\overset{―}{w_{h}}}^{E}$ . If $E$ is an inner edge, then the mean value on the corresponding edge of the neighbor element is the same. The normal derivative $\frac{\partial u}{\partial n}$ on the neighbor element is (up to the sign) the same. If $E$ is an edge on the Dirichlet boundary, then the mean value is 0. This allows to subtract edge mean values:

r (w_{h}) = \sum_{T} \sum_{E \subset T} \int_{E} \frac{\partial u}{\partial n} (w_{h} - {\overset{―}{w_{h}}}^{E}) d s

Since $\int_{E} w_{h} - {\overset{―}{w_{h}}}^{E} d s = 0$ , we may insert the constant function $\frac{\partial I_{h} u}{\partial n}$ on each edge:

r (w_{h}) = \sum_{T} \sum_{E \subset T} \int_{E} (\frac{\partial u}{\partial n} - \frac{\partial I_{h} u}{\partial n}) (w_{h} - {\overset{―}{w_{h}}}^{E}) d s

Apply Cauchy-Schwarz on $L_{2} (E)$ :

r (w_{h}) = \sum_{T} \sum_{E \subset T} ‖ \nabla (u - I_{h} u) ‖_{L_{2} (E)} ‖ w_{h} - {\overset{―}{w_{h}}}^{E} ‖_{L_{2} (E)}

To estimate these terms, we transform to the reference element $\hat{T}$ , where we apply the Bramble Hilbert lemma. Let $T = Φ_{T} (\hat{T})$ , and set

\hat{u} = u \circ Φ_{T} {\hat{w}}_{h} = w_{h} \circ Φ_{T}

There hold the scaling estimates

\begin{array}{rcl} | w_{h} |_{H^{1} (T)} & ≃ & | {\hat{w}}_{h} |_{H^{1} (\hat{T})} \\ ‖ w_{h} - {\overset{―}{w_{h}}}^{E} ‖_{L_{2} (E)} & ≃ & h_{E}^{1 / 2} ‖ {\hat{w}}_{h} - {\overset{―}{{\hat{w}}_{h}}}^{\hat{E}} ‖_{L_{2} (\hat{E})} \\ | u |_{H^{2} (T)} & ≃ & h_{T}^{- 1} | \hat{u} |_{H^{2} (\hat{T})} \\ ‖ \nabla (u - I_{h} u) ‖_{L_{2} (E)} & ≃ & h_{E}^{- 1 / 2} ‖ \nabla (\hat{u} - {\hat{I}}_{h} \hat{u}) ‖_{L_{2} (E)} \end{array}

On the reference element, we apply the Bramble Hilbert lemma, once for $w_{h}$ , and once for $u$ . The linear operator

L : H^{1} (\hat{T}) \to L_{2} (\hat{E}) : {\hat{w}}_{h} \to {\hat{w}}_{h} - {\overset{―}{{\hat{w}}_{h}}}^{\hat{E}}

is bounded on $H^{1} (\hat{T})$ (trace theorem), and $L w = 0$ for $w \in P_{0}$ , thus

‖ {\hat{w}}_{h} - {\overset{―}{{\hat{w}}_{h}}}^{\hat{E}} ‖_{L_{2} (\hat{E})} ⪯ | {\hat{w}}_{h} |_{H^{1} (\hat{T})}

Similar for the term in $u$ : There is $‖ \nabla (u - I_{h} u) ‖_{L_{2} (\hat{E})} ⪯ ‖ u ‖_{H^{2} (\hat{T})}$ , and $u - I_{h} u$ vanishes for $u \in P^{1}$ .

Rescaling to the element $T$ leads to

\begin{array}{rcl} ‖ w_{h} - {\overset{―}{w_{h}}}^{E} ‖_{L_{2} (E)} & ⪯ & h^{1 / 2} | w_{h} |_{H^{1} (T)} \\ ‖ \nabla (u - I_{h} u) ‖_{L_{2} (E)} & ⪯ & h^{1 / 2} | u |_{H^{2} (T)} \end{array}

This bounds the consistency term

r (w_{h}) ⪯ \sum_{T} h | u |_{H^{2} (T)} | w_{h} |_{H^{1} (T)} ⪯ h ‖ u ‖_{H^{2} (Ω)} ‖ w_{h} ‖_{h} .

Thus, the second lemma of Strang provides the error estimate

‖ u - u_{h} ‖ ⪯ h ‖ u ‖_{H^{2}}

There are several applications where the non-conforming $P^{1}$ triangle is of advantage:

The $L_{2}$ matrix is diagonal (exercises)
It can be used for the approximation of problems in fluid dynamics described by the Navier Stokes equations (see later).
The finite element matrix has exactly 5 non-zero entries in each row associated with inner edges. That allows simplifications in the matrix implementation.

25.3. Solving Stokes’ equation with the non-conforming $P^{1}$ -triangle#

The weak formulation of Stokes` equation is: Find the velocity vector-field $u \in [H^{1}]^{2}$ and a pressure $p \in L_{2}$ such that

\begin{array}{r} \begin{array}{ccccll} (\nabla u, \nabla v) & + & (div v, p) & = & (f, v) & \forall v \in [H^{1}]^{2} \\ (div u, q) & = & 0 & \forall q \in L_{2} \end{array} \end{array}

The first equation is equilibrium of forces, $- Δ u = f - \nabla p$ , the second one is the incompressibility constraint $div u = 0$ .

We will discuss variational formulation with constraints later in the chapter Mixed finite element methods. A key property is the following

Lemma: There holds the commuting diagram property:

$div I_{T} u = Π^{P^{0}} div u$

Here, $I_{T}$ is the interpolation operator by mean-values on edges into the non-conforming $P^{1}$ space, and $Π^{P^{0}}$ is the $L_{2}$ -orthogonal projection onto constants.

Proof: Both sides are constant functions on each triangle $T$ . We apply Gauss` theorem to show that they are the same:

\begin{aligned} \int_{T} div I_{T} u & = \int_{\partial T} (I_{T} u)_{n} = \sum_{E \in \partial T} \int_{E} (I_{T} u)_{n} \\ = \sum_{E \in \partial T} \int_{E} u_{n} = \int_{\partial T} u_{n} = \int_{T} div u \\ = \int_{T} Π^{P^{0}} div u \end{aligned}

from ngsolve import *
from netgen.occ import *
from ngsolve.webgui import Draw

rect = Rectangle(10,1).Face()
rect.edges.name = "wall"
rect.edges.Min(X).name = "inlet"
rect.edges.Max(X).name = "outlet"

mesh = Mesh(OCCGeometry(rect, dim=2).GenerateMesh(maxh=0.2))
Draw (mesh);

Build a product spaces with components $\vec{u}$ and $p$ . Express $div u = tr \nabla u$ :

# fesh1 = H1(mesh, order=2, dirichlet="wall|inlet")
fesh1 = FESpace("nonconforming", mesh, order=1, dirichlet="wall|inlet")
fesl2 = L2(mesh, order=0)
fes = VectorValued(fesh1) * fesl2

u,p = fes.TrialFunction()
v,q = fes.TestFunction()

a = BilinearForm(fes)
a += InnerProduct(grad(u),grad(v))*dx
a += Trace(grad(u))*q*dx + Trace(grad(v))*p*dx 
a.Assemble()

gfu = GridFunction(fes)
gfu.components[0].Set((y*(1-y), 0), definedon=mesh.Boundaries("inlet"))

gfu.vec.data -= a.mat.Inverse(fes.FreeDofs())@a.mat * gfu.vec

Draw (gfu.components[0], mesh, vectors={"grid_size":80})
Draw (gfu.components[1]);

Non-conforming Finite Element Methods

Contents

25. Non-conforming Finite Element Methods#

25.1. The First Lemma of Strang#

25.2. The Second Lemma of Strang#

25.3. Solving Stokes’ equation with the non-conforming P1-triangle#

25.3. Solving Stokes’ equation with the non-conforming $P^{1}$ -triangle#