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

\[ \mbox{Find } u_h \in V_h: \qquad A(u_h, v_h) = f(v_h) \qquad \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 \rightarrow {\mathbb R} \]

and

\[ f_h(.) : V_h \rightarrow {\mathbb 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 \(\alpha_1\) independent of the mesh-size such that

\[ A_h (v_h, v_h) \geq \alpha_1 \, \| v_h \|_V^2 \qquad \forall \, v_h \in V_h \]

The finite element problem is defined as

\[ \mbox{Find } u_h \in V_h: \qquad A_h (u_h, v_h) = f_h (v_h) \qquad \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{eqnarray*} \| u - u_h \| & \preceq & \inf_{v_h \in V_h} \left\{ \| u - v_h \| + \sup_{w_h \in V_h} \frac{|A(v_h, w_h) - A_h (v_h, w_h)|}{\| w_h \|} \right\} \\ & & + \sup_{w_h \in V_h} \frac{f(w_h) - f_h (w_h)}{\| w_h \|} \end{eqnarray*}\]

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{align*} \alpha_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{align*}\]

Divide by \(\| u_h - v_h \| = \| w_h \|\), and use the continuity of \(A(.,.)\):

\[ \| u_h - v_h \| \preceq \| 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{align*} \| u - u_h \| & \leq \inf_{v_h \in V_h} \| u - v_h \| + \| v_h - u_h \| \\ & \preceq \inf_{v_h \in V_h} \left\{ \| 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 \| } \right\} \\ & \preceq \inf_{v_h \in V_h} \left\{ \| 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 \| } \right\} \end{align*}\]

The lemma is proven. \(\Box\)

Example: Lumping of the \(L_2\) bilinear-form: Define the \(H^1\) - bilinear-form

\[ A(u,v) = \int_\Omega \nabla u \cdot \nabla v + \int_\Omega u v \, dx, \]

and perform Galerkin discretization with \(P^1\) triangles. The second term leads to a non-diagonal matrix. The vertex integration rule

\[ \int_T v \, dx \approx \frac{|T|}{3} \sum_{\alpha = 1}^3 v(x_{T,\alpha}) \]

is exact for \(v \in P^1\). We apply this integration rule for the term \(\int u v \, dx\):

\[ A_h(u,v) = \int \nabla u \cdot \nabla v + \sum_{T \in {\cal T}} \frac{|T|}{3} \sum_{\alpha = 1}^3 u(x_{T,\alpha}) v(x_{T,\alpha}) \]

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 \(\varphi_i\), we obtain a diagonal matrix for the second term:

\[\begin{split} \varphi_i (x_{T,\alpha}) \varphi_j (x_{T,\alpha}) = \left\{ \begin{array}{cl} 1 & \mbox{for } x_i = x_j = x_{T,\alpha} \\ 0 & \mbox{else} \end{array} \right. \end{split}\]

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

\[ \sum_T \frac{|T|}{3} \sum_{\alpha = 1}^3 |v_h(x_{T,\alpha})|^2 \geq \alpha_1 \sum_T \int_T | v_h |^2 \, dx \qquad \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 \, dx - \frac{|T|}{3} \sum_{\alpha=1}^3 u_h(x_\alpha) v_h(x_\alpha) | \preceq 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) \preceq h^2 \| u_h \|_{H^1(\Omega)} \, \| v_h \|_{H^1(\Omega)} \]

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

\[\begin{align*} \| u - u_h \|_{H^1} & \preceq \inf_{v_h \in V_h} \left\{ \| 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}} \right\} \\ & \preceq \| 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}} \\ & \preceq 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}} \\ & \preceq h \, \| u \|_{H^2} \end{align*}\]

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 \alpha_1 \| v_h \|_h^2 \qquad \forall \, v_h \in V_h \]

• continuity:

  \[ A_h (u, v_h) \leq \alpha_2 \| u \|_h \| v_h \|_h \qquad \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 \preceq \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{align*} \alpha_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{align*}\]

Again, divide by \(\| u_h - v_h\|\), and use continuity of \(A_h(.,.)\):

\[ \| u_h - v_h \|_h \preceq \| 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. \(\Box\)

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^{nc} := \{ v \in L_2 : v_{|T} \in P^1(T), \mbox{and }v \mbox{ is continuous in edge mid-points} \} \]

The functions in \(V_h^{nc}\) are not continuous across edges, and thus, \(V_h^{nc}\) 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 {\cal T}} \int_T \nabla u \nabla v \, dx \qquad \forall \, u, v \in V + V_h^{nc} \]

and

\[ \| v \|_h^2 := \sum_{T \in {\cal T}} \| \nabla v \|_{L_2(T)}^2 \qquad \forall \, v \in V + V_h^{nc} \]

We consider the Dirichlet-problem with \(u = 0\) on \(\Gamma_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^{nc}\). Let \(I_h : H^2 \rightarrow 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^{nc}\):

\[ \inf_{v_h \in V_h^{nc}} \| u - v_h \|_h \leq \| u - I_h u \|_{H^1} \preceq h \, \| u \|_{H^2} \]

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

\[\begin{eqnarray*} 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 \, dx \\ & = & \sum_T \int_{\partial T} \frac{\partial u}{\partial n} w_h \, ds - \sum_T \int_T (\Delta u + f) \, w_h \, ds \\ & = & \sum_T \int_{\partial T} \frac{\partial u}{\partial n} w_h \, ds \end{eqnarray*}\]

Let \(E\) be an edge of the triangle \(T\). Define the mean value \(\overline{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 - \overline{w_h}^E) \, ds \]

Since \(\int_E w_h - \overline{w_h}^E \, ds = 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 \left(\frac{\partial u}{\partial n} - \frac{\partial I_h u}{\partial n} \right) (w_h - \overline{w_h}^E) \, ds \]

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 - \overline{w_h}^E \|_{L_2(E)} \]

To estimate these terms, we transform to the reference element \(\widehat T\), where we apply the Bramble Hilbert lemma. Let \(T = \Phi_T(\widehat T)\), and set

\[ \widehat u = u \circ \Phi_T \qquad \widehat w_h = w_h \circ \Phi_T \]

There hold the scaling estimates

\[\begin{eqnarray*} | w_h |_{H^1(T)} & \simeq & | \widehat w_h |_{H^1(\widehat T)} \\ \| w_h - \overline {w_h}^E \|_{L_2(E)} & \simeq & h_E^{1/2} \| \widehat w_h - \overline{\widehat w_h}^{\widehat E} \|_{L_2(\widehat E)} \\ | u |_{H^2(T)} & \simeq & h_T^{-1} | \widehat u |_{H^2(\widehat T)} \\ \| \nabla (u - I_h u) \|_{L_2(E)} & \simeq & h_E^{-1/2} \| \nabla (\widehat u - \widehat I_h \widehat u) \|_{L_2(E)} \end{eqnarray*}\]

On the reference element, we apply the Bramble Hilbert lemma, once for \(w_h\), and once for \(u\). The linear operator

\[ L : H^1(\widehat T) \rightarrow L_2(\widehat E) : \widehat w_h \rightarrow \widehat w_h - \overline {\widehat w_h}^{\widehat E} \]

is bounded on \(H^1(\widehat T)\) (trace theorem), and \(L w = 0\) for \(w \in P_0\), thus

\[ \| \widehat w_h - \overline{\widehat w_h}^{\widehat E} \|_{L_2(\widehat E)} \preceq | \widehat w_h |_{H^1(\widehat T)} \]

Similar for the term in \(u\): There is \(\| \nabla (u - I_h u) \|_{L_2(\widehat E)} \preceq \| u \|_{H^2(\widehat T)}\), and \(u - I_h u\) vanishes for \(u \in P^1\).

Rescaling to the element \(T\) leads to

\[\begin{eqnarray*} \| w_h - \overline{w_h}^E \|_{L_2(E)} & \preceq & h^{1/2} \, | w_h |_{H^1(T)} \\ \| \nabla (u - I_h u) \|_{L_2(E)} & \preceq & h^{1/2} \, | u |_{H^2(T)} \end{eqnarray*}\]

This bounds the consistency term

\[ r(w_h) \preceq \sum_T h \, | u |_{H^2(T)} | w_h |_{H^1(T)} \preceq h \, \| u \|_{H^2(\Omega)} \, \| w_h \|_h. \]

Thus, the second lemma of Strang provides the error estimate

\[ \| u - u_h \| \preceq 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{split} \begin{array}{ccccll} (\nabla u, \nabla v) & + & (\operatorname{div} v, p) & =& (f,v) & \forall \, v \in [H^1]^2 \\ (\operatorname{div} u,q) & & & = & 0 & \forall \, q \in L_2 \end{array} \end{split}\]

The first equation is equilibrium of forces, \(-\Delta u = f - \nabla p\), the second one is the incompressibility constraint \(\operatorname{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:

\[ \operatorname{div} I_T u = \Pi ^{P^0} \operatorname{div} u \]

Here, \(I_T\) is the interpolation operator by mean-values on edges into the non-conforming \(P^1\) space, and \(\Pi ^{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{align*} \int_T \operatorname{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 \operatorname{div} u \\ & = \int_T \Pi ^{P^0} \operatorname{div} u \end{align*}\]

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 \(\operatorname{div} u = \operatorname{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]);