Crouzeix’s divergence free basis for the Navier-Stokes equation

July 28, 2021 15 minute read

This is the solution (long read) to a Chapter 11 exercise of the book by Brenner and Scott (Brenner and Scott, 2008). Personally I really like this problem, as it trains us to think about the orientation of edges locally in an element versus the global orientation in a mesh.

Recently, I am reading some constructions related to a divergence-free basis for Stokes flow using virtual element, in a sense, it resembles a lot with this problem’s “macro-element” style approach. So, here I revisited my proof 10 years ago and post it here.

An interesting problem on verifying a divergence free basisPermalink

Let

$\boldsymbol{V}_h^*$ be defined as follows:

$\begin{aligned} \boldsymbol{V}_h^* = & \{\boldsymbol{v}\in \boldsymbol{L}^2(\Omega): v \text{ is linear for all } K\in \mathcal{T}_h, \\ & v \text{ is continuous at the midpoints of inter-element boundaries,} \\ & v \text{ is zero on any edge in } \partial \Omega \}. \end{aligned}$ Let

$\boldsymbol{W}_h := \{ \boldsymbol{v} \in \boldsymbol{V}_h^*: \mathrm{div}_h \boldsymbol{v} = 0 \}$ Denote by

$\psi_e$ the piecewise linear function that takes the value

$1$ at the midpoint of the edge

$e$ and

$0$ at all other midpoints. The first kind of basis function of

$\boldsymbol{W}_h$ are associated with internal edges. Let

$e$ be an internal edge of

$\mathcal{T}_h$ and

$\boldsymbol{t}_e$ be a unit vector tangential to

$e$ , then

$\boldsymbol{\phi}_e := \psi_e \boldsymbol{t}_e$ . The second kind of basis functions are associated with internal vertices. Let

$p$ be an internal vertex and let

$e_1, \dots, e_l$ be the edge in

$\mathcal{T}_h$ that have

$p$ as an endpoint, then

$\boldsymbol{\phi}_p := \sum^l_{i=1}\frac{1}{|e_i|} \psi_{e_i} \boldsymbol{n}_{e_i}$ , where

$\boldsymbol{n}_{e_i}$ is the unit normal vector to

$e_i$ , pointing in the counterclockwise direction with respect to

$p$ (see the figure below).

divergence free basis flow

Show that a basis for

$\boldsymbol{W}_h$ is given by the union of the two sets:

$\{\boldsymbol{\phi}_e: e \text{ is an internal edge of } \mathcal{T}_h\}$ and

$\{\boldsymbol{\phi}_p: e \text{ is an internal vertex of } \mathcal{T}_h\}.$

ProofsPermalink

Part 1: local divergence-freeness for basisPermalink

For an edge basis $\boldsymbol{\phi}_e = \psi_e \boldsymbol{t}_e$ , consider any $T\subset \omega_e$

$\operatorname{div} \boldsymbol{\phi}_e\big|_T = \operatorname{div}(\psi_e \boldsymbol{t}_e) = \nabla \psi_e \cdot \boldsymbol{t}_e = \boldsymbol{n}_e \cdot \boldsymbol{t}_e \frac{|\tilde{e}|}{2|\tilde{T}|} = 0$

where $\tilde{T}$ is the small triangle connecting three midpoints in $T$ .

Now consider the vertex basis $\boldsymbol{\phi}_p$ on some element $T\in\omega_p$ , denote the joint edge of $p$ as $e_1$ and $e_2$ we have

$\operatorname{div} \boldsymbol{\phi}_p\big\vert_T = \frac{1}{|e_1|} \nabla \psi_{e_1} \cdot \boldsymbol{n}_{e_1} + \frac{1}{|e_2|} \nabla \psi_{e_2} \cdot \boldsymbol{n}_{e_2}.$

Refer to the figure above, $\nabla \psi_{e_1}$ is the inner normal vector of $\tilde{e}_1$ with respect to the triangle $\tilde{T}$ , denote the inner normal vector as $\boldsymbol{\nu}_{\tilde{e}_1}$ then we have

$\nabla \psi_{e_1} = \frac{1}{2|\tilde{T}|} \boldsymbol{\nu}_{\tilde{e}_1} = -\frac{1}{2|\tilde{T}|} |\tilde{e}_1| \boldsymbol{n}_{e_1} = -\frac{|e_1|}{|T|}\boldsymbol{n}_{e_1},$

therefore

$\frac{1}{|e_1|} \nabla \psi_{e_1} \cdot \boldsymbol{n}_{e_1} = -\frac{1}{|T|}.$

Similarly, we have

$\frac{1}{|e_2|} \nabla \psi_{e_2} \cdot \boldsymbol{n}_{e_2} = \frac{1}{|T|},$

hence

$\operatorname{div} \boldsymbol{\phi}_p\big|_T = 0.$

Part 2: write a vector as the linear combination of basis functionsPermalink

Now for $\forall M$ which is the midpoint of an edge $e_i$ in this patch, consider $T^+$ and $T^-$ , where $\omega_e = T^+ \cup T^-$ . For the space $\boldsymbol{V}_h^*$ in which the functions are continuous at $M$ , consider $\forall \boldsymbol{w} \in \boldsymbol{V}_h^*$ , $\boldsymbol{w}(M)$ is a constant vector, hence we could decompose $\boldsymbol{w}$ into two orthogonal vectors in the directions of $\boldsymbol{t}_e$ and $\boldsymbol{n}_e$ , respectively. At point $M$ we have:

$\boldsymbol{w}(M) = (\boldsymbol{w}\cdot\boldsymbol{t}_e)\boldsymbol{t}_e + (\boldsymbol{w}\cdot\boldsymbol{n}_e)\boldsymbol{n}_e.$

This is true because at point $M$ we could express $\boldsymbol{w} = |\boldsymbol{w}|\boldsymbol{\xi}$ , where $\boldsymbol{\xi}$ is the unit vector in $\boldsymbol{w}$ ’s direction, denote the angle between $\boldsymbol{w}$ and $\boldsymbol{t}_e$ as $\theta$ ,

$\boldsymbol{\xi} = \boldsymbol{t}_e\cos \theta + \boldsymbol{n}_e\sin \theta= \boldsymbol{t}_e \left(\frac{\boldsymbol{w}}{|\boldsymbol{w}|}\cdot\boldsymbol{t}_e\right) + \boldsymbol{n}_e\left(\frac{\boldsymbol{w}}{|\boldsymbol{w}|}\cdot\boldsymbol{n}_e\right).$

It is easy to see that $\boldsymbol{\xi}\cdot \boldsymbol{\xi} = 1$ since $\boldsymbol{t}_e$ and $\boldsymbol{n}_e$ are orthogonal. Now at $M$ we know $\psi_e = 1$ , hence we could express at $M$ :

$\boldsymbol{w} = (\boldsymbol{w}\cdot\boldsymbol{t}_e)\psi_e\boldsymbol{t}_e + (\boldsymbol{w}\cdot\boldsymbol{n}_e)\psi_e\boldsymbol{n}_e.$

Since for all other four midpoints within the $\omega_e = T^+ \cup T^-$ , the basis functions $\psi_{e’}$ for $e’\subset \partial \omega_e$ are all zero at $M$ , and linear within the triangle, any linear vector that could be decomposed this way at $M$ bears the same continuity with two basis function $\psi_e \boldsymbol{t}_e$ and $\psi_e \boldsymbol{n}_e$ at the same edge. Other edge bases vanish at $M$ . In other words, the linear vector could be written in this way is continuous at $M$ .

Now locally within any element $T$ that has edge $e_1$ , $e_2$ , and $e_3$ , with midpoints $M_1$ , $M_2$ and $M_3$ , $\forall x\in T$ , for any linear vector $\boldsymbol{w}$ decompose

$\begin{aligned} \boldsymbol{w}(x) &= \boldsymbol{w}(\psi_{e_1}(x)M_1 + \psi_{e_2}(x)M_2 + \psi_{e_3}(x)M_3 ) \\ &= \psi_{e_1}(x) \boldsymbol{w}(M_1) + \psi_{e_2}(x) \boldsymbol{w}(M_2) + \psi_{e_3}(x) \boldsymbol{w}(M_3) \end{aligned}$

for each $\psi_{e_j}(x)\boldsymbol{w}(M_j)$ which is a constant at $x$ , decompose it in the same fashion as before along the direction of $\boldsymbol{t}_{e_j}$ and $\boldsymbol{n}_{e_j}$ . Hence $\forall \boldsymbol{w}\in \boldsymbol{V}_h^*$ , globally we could write it as

$\boldsymbol{w} = \sum_{e\in \mathcal{E}_h} \Bigl((\boldsymbol{w}(M_e)\cdot \boldsymbol{t}_e) \boldsymbol{\phi}_e + (\boldsymbol{w}(M_e)\cdot \boldsymbol{n}_e) \boldsymbol{\eta}_e \Bigr)$

where $\boldsymbol{\phi}_e = \psi_e \boldsymbol{t}_e$ and $\boldsymbol{\eta}_e = \psi_e \boldsymbol{n}_e$ . Moreover, $\boldsymbol{\phi}_e\perp \boldsymbol{\eta}_e$ , i.e., they are linearly independent. Now we have shown $\{\boldsymbol{\phi}_e\}$ and $\{\boldsymbol{\eta}_e\}$ form a set of basis functions for $\boldsymbol{V}^*_h$ .

Part 3: more on the divergence-freenessPermalink

For the divergence zero subspace, we have already shown

$\operatorname{div}_h \boldsymbol{\phi}_e = 0.$

It suffices to show that $\forall \boldsymbol{v}$ such that

$\boldsymbol{v} = \sum_{e\in \mathcal{E}_h} (\boldsymbol{w}(M_e)\cdot \boldsymbol{n}_e) \boldsymbol{\eta}_e = \sum_{e\in \mathcal{E}_h} b_e \boldsymbol{\eta}_e\in \boldsymbol{V}^*_h$

of which $\operatorname{div}_h\boldsymbol{v} = 0$ . This $\boldsymbol{v}$ can be written as

$\boldsymbol{v} = \sum_{p\in \mathcal{V}_h} b_p \boldsymbol{\phi}_p.$

Now consider within an arbitrary element $T\subset \omega_e$ , and assign a fixed unit outer normal vector to each edge of $T$ as in the following figure we have:

$\begin{aligned} 0 = \operatorname{div}_h\boldsymbol{v}\Big|_T &= \sum_{e\in \mathcal{E}_h} \Bigl( (\boldsymbol{v}(M_e)\cdot \boldsymbol{n}_e) \operatorname{div}_h(\boldsymbol{\eta}_e) + \nabla_h (\boldsymbol{v}(M_e)\cdot \boldsymbol{n}_e) \boldsymbol{\eta}_e \Bigr) \\ 0&= \sum_{e\subset \partial T} (\boldsymbol{v}(M_e)\cdot \boldsymbol{n}_e) (\nabla \psi_e \cdot \boldsymbol{n}_e) \quad \\ 0&= (\boldsymbol{v}(M_1)\cdot \boldsymbol{n}_{e_1}) \frac{|e_1|}{|T|} + (\boldsymbol{v}(M_2)\cdot \boldsymbol{n}_{e_2}) \frac{|e_2|}{|T|} + (\boldsymbol{v}(M_3)\cdot \boldsymbol{n}_{e_3}) \frac{|e_3|}{|T|} \end{aligned}$

rewriting above using the coefficients we have

$b_{e_1}|e_1| + b_{e_2}|e_2| + b_{e_3}|e_3| = 0. \tag{1}\label{divfree}$

Therefore, as long as the coefficients of edge basis functions satisfy above relation, the vector field is divergence free within this element $T$ .

Now consider the linear vector field $\boldsymbol{u} = \sum_{p\in \mathcal{V}_h} b_p \boldsymbol{\phi}_p$ within $T$ , in the following equation we would like to rewrite $\boldsymbol{\phi}_p$ ’s counterclockwise normal vector aligning with the local unit outward normal, and we have:

$\begin{aligned} \boldsymbol{u} &= b_{p_1} \boldsymbol{\phi}_{p_1} + b_{p_2} \boldsymbol{\phi}_{p_2} + b_{p_3} \boldsymbol{\phi}_{p_3} \\ &= b_{p_1} \Bigl(\frac{1}{|e_3|}\boldsymbol{\eta}_{e_3} - \frac{1}{|e_2|}\boldsymbol{\eta}_{e_2} \Bigr) + b_{p_2} \Bigl(\frac{1}{|e_1|}\boldsymbol{\eta}_{e_1} - \frac{1}{|e_3|}\boldsymbol{\eta}_{e_3}\Bigr) + b_{p_3} \Bigl(\frac{1}{|e_2|}\boldsymbol{\eta}_{e_2} - \frac{1}{|e_1|}\boldsymbol{\eta}_{e_1} \Bigr) \\ &= \frac{1}{|e_1|}(b_{p_2} - b_{p_3})\boldsymbol{\eta}_{e_1} + \frac{1}{|e_2|}(b_{p_3} - b_{p_1})\boldsymbol{\eta}_{e_2} + \frac{1}{|e_3|}(b_{p_1} - b_{p_2})\boldsymbol{\eta}_{e_3} \\ &:= b_{e_1}\boldsymbol{\eta}_{e_1} + b_{e_2}\boldsymbol{\eta}_{e_2} + b_{e_3}\boldsymbol{\eta}_{e_3} \end{aligned}$

It is obvious that here the coefficients satisfy $\eqref{divfree}$ . Also the converse is also true, for a divergence free vector

$\boldsymbol{u} = b_{e_1}\boldsymbol{\eta}_{e_1} + b_{e_2}\boldsymbol{\eta}_{e_2} + b_{e_3}\boldsymbol{\eta}_{e_3}$

within element $T$ , by $\eqref{divfree}$ we have

$\begin{aligned} \boldsymbol{u} &= b_{e_1}\boldsymbol{\eta}_{e_1} + b_{e_2}\boldsymbol{\eta}_{e_2} - \frac{1}{|e_3|} (b_{e_1}|e_1| + b_{e_2}|e_2|)\boldsymbol{\eta}_{e_3} \\ &= b_{e_1}|e_1|\bigl(\frac{1}{|e_1|}\boldsymbol{\eta}_{e_1} - \frac{1}{|e_3|}\boldsymbol{\eta}_{e_3} \bigr) + b_{e_2}|e_2|\bigl(\frac{1}{|e_2|}\boldsymbol{\eta}_{e_2} - \frac{1}{|e_3|}\boldsymbol{\eta}_{e_3} \bigr) \\ &= b_{p_2} \boldsymbol{\phi}_{p_2} + b_{p_1}\boldsymbol{\phi}_{p_1} \end{aligned}$

or a linear combination of $\boldsymbol{\phi}_{p_3}$ and $\boldsymbol{\phi}_{p_2}$ , or $\boldsymbol{\phi}_{p_3}$ and $\boldsymbol{\phi}_{p_1}$ due to relation $\eqref{divfree}$ . Hence here we proved that locally within each element $T$ , any divergence free linear vector field could be written as the linear combination of the $\{\boldsymbol{\phi}_e\}_{e\subset \partial T}$ and $\{\boldsymbol{\phi}_p\}_{p\in \mathcal{V}_T}$ .

Part 4: how does the global continuity condition work?Permalink

Now we would like to prove globally glued together with the continuity condition, $\boldsymbol{\phi}_p$ will bear the following form globally:

$\boldsymbol{\phi}_p = \sum^l_{i=1} \frac{1}{|e_i|} \psi_{e_i} \boldsymbol{n}_{e_i}$

Consider the vertex $p$ , the edge $e_1$ , $e_2$ , and $e_3$ , two elements $T_1$ and $T_2$ , unit normal for each edge fixed as in the following figure, and $\boldsymbol{\eta}_{e_i} = \psi_{e_i}\boldsymbol{n}_{e_i}$ :

We already knew that locally in $T_1$ :

$\boldsymbol{\phi}_p\big|_{T_1} = \frac{1}{|e_1|}\, \boldsymbol{\eta}_{e_1}\big|_{T_1} + \frac{1}{|e_2|} \,\boldsymbol{\eta}_{e_2}\big|_{T_1}$

clearly coefficients satisfy $\eqref{divfree}$ . Similarly, we could get

$\boldsymbol{\phi}_p\big|_{T_2} = \frac{1}{|e_2|} \,\boldsymbol{\eta}_{e_2}\big|_{T_2} + \frac{1}{|e_3|} \, \boldsymbol{\eta}_{e_3}\big|_{T_2}$

Consider any $\boldsymbol{u}$ that is divergence free within $T_1$ and $T_2$ , that

$\boldsymbol{u}\big|_{T_1} = a_1 \boldsymbol{\phi}_p\big|_{T_1}, \quad \boldsymbol{u}\big|_{T_2} = a_2 \boldsymbol{\phi}_p\big|_{T_2}$

and while continuous at $M$ :

$\lim_{x\to M,x\in T_1} \boldsymbol{u}= \frac{a_1}{|e_2|} \boldsymbol{n}_{e_2} = \lim_{x\to M,x\in T_2} \boldsymbol{u} = \frac{a_2}{|e_2|} \boldsymbol{n}_{e_2}$

this implies $a_1 = a_2$ , i.e., the coefficient for each $\boldsymbol{\phi}_{p}\big|_{T_i}$ should be the same when marching counterclockwisely around this patch, where $T_i \subset \omega_p$ is an element within this patch. Together with last part, this implies

$\{\boldsymbol{\phi}_p\} = \{\boldsymbol{v}: \boldsymbol{v} = \sum_{e\in \mathcal{T}_h} b_e\boldsymbol{\eta}_e\}\bigcap \{\boldsymbol{v}: \operatorname{div}_h \boldsymbol{v} = 0\}$

Therefore, $\{\boldsymbol{\phi}_e\}$ and $\{\boldsymbol{\phi}_p\}$ are shown to be the basis functions for the space

$\boldsymbol{W}_h = \{\boldsymbol{v}\in \boldsymbol{V}_h^*: \operatorname{div}_h \boldsymbol{v} = 0\}$ .

A final remarkPermalink

Notice for general vector field $\boldsymbol{u}$ , the nodal interpolation of it onto $\boldsymbol{V}^*_h$ is

$\operatorname{I}_{\boldsymbol{V}^*_h}\boldsymbol{u} = \sum_{T\in \mathcal{T}_h}\sum_{e\subset \partial T} \Bigl((\boldsymbol{u},\boldsymbol{\phi}_e)_T \boldsymbol{\phi}_e + (\boldsymbol{u}, \boldsymbol{\eta}_e)_T \boldsymbol{\eta}_e \Bigr).$

and this projection does not coincide with the “nodal” interpolation of linear vector field $\boldsymbol{u}$ using pointwise value as coefficient, since

$\begin{aligned} (\boldsymbol{u},\boldsymbol{\phi}_e)_T &= \int_T \boldsymbol{u} \cdot \boldsymbol{t}_e \psi_e \,dx \\ &= \frac{|e|}{2|T|} \int_T \boldsymbol{u} \cdot \nabla^{\perp} \psi_e^2 \,dx \\ &= \frac{|e|}{2|T|} \int_T \nabla\times \boldsymbol{u} \,\psi_e^2\,dx + \frac{|e|}{2|T|} \int_{\partial K} \boldsymbol{u}\cdot \boldsymbol{t} \psi_e^2 \,ds. \end{aligned}$

ReferencesPermalink

Brenner, S. C. and Scott, R. (2008) The mathematical theory of finite element methods. Springer.
```
@book{BrennerScott2008mathematical,
  author = {Brenner, Susanne C and Scott, Ridgway},
  title = {The mathematical theory of finite element methods},
  publisher = {Springer},
  volume = {15},
  date-modified = {2017-11-02 17:52:56 +0000},
  year = {2008}
}
```

Share on

X Facebook LinkedIn Bluesky

Crouzeix’s divergence free basis for the Navier-Stokes equation

An interesting problem on verifying a divergence free basisPermalink

ProofsPermalink

Part 1: local divergence-freeness for basisPermalink

Part 2: write a vector as the linear combination of basis functionsPermalink

Part 3: more on the divergence-freenessPermalink

Part 4: how does the global continuity condition work?Permalink

A final remarkPermalink

ReferencesPermalink

Share on

Comments

You May Also Enjoy

Neural Finite Element part II: Einstein summation is all we need

Neural Finite Element part I: how to implement MATLAB’s accumarray in PyTorch

Amateur mathematician studies Quantum Field Theory part 2: Maxwell equations

Amateur mathematician studies Quantum Field Theory part I