# Why aren't face integrals for an element calculated in FEM but they show up in FVM?

Consider the Laplace problem: \begin{align} -\nabla^2 u = f \qquad \text{in } \Omega \\ u = 0 \qquad \text{on } \Gamma \end{align} The weak problem is find $$u_h \in V \subset H^1$$ such that $$\forall v_h \in K \subset H^1_0$$, $$u_h$$ satisfies:
$$$$\int_\Omega \nabla v_h \cdot \nabla u_h \, d\Omega = \int_\Omega f v_h \, d\Omega$$$$ FVM is said to be a special case of FEM where the space of test functions is $$1$$. In FVM we get to work with the integral $$$$\oint_e \nabla u \, d\Omega$$$$ for each local element. Why does this integral show up in FVM but not in Galerkin FEM? This doesn't feel consistent to me.

The solution $$u\in H^1(\Omega)$$ for the Poisson problem is continuous-differentiable and therefore we can apply the div-theorem such that we end up with the weak form

\begin{align} -\int_\Omega \Delta u v dx= \int_\Omega \nabla u \nabla v \ dx - \underbrace{\int_{\Gamma} \nabla u \cdot n v \ ds}_{=0\text{ since } u=0 \text{ on } \Gamma} =\int_\Omega f v \ dx \quad \forall v\in H^1_0(\Omega) \end{align}

As the comments suggest you need here that the weak-gradients $$\nabla u, \nabla v$$ exist, which is why we at least need that $$u,v\in H^1(\Omega)$$.

For a triangulation $$\mathcal{T}$$ with elements $$T\in \mathcal{T}$$ you note that

\begin{align} \sum_{T\in \mathcal{T}} \int_{\partial T}\nabla u \cdot n v \ ds=\int_{\Gamma} \nabla u\cdot n v \ ds = 0 \end{align} must be true. Therefore all the contribution on the element-faces inside the domain $$\Omega$$ musst cancel out, or to be precise, the normal fluxes $$\nabla u \cdot n$$ are continuous at the element boundaries.

Now we approximate the solution of the weak-form above by using discrete subspaces $$V_h\subset H^1(\Omega)$$, $$V_{h0}\subset H^1_0(\Omega)$$. If you rewrite the weak form above for the discrete setting we therefore get:

\begin{align} \int_\Omega \nabla u_h \nabla v_h \ dx - \underbrace{\int_{\Gamma} \nabla u_h \cdot n v_h \ ds}_{=0\text{ since } u_h=0 \text{ on } \Gamma}\quad =\int_\Omega f v_h \ dx \quad \forall v_h\in V_{h0}(\Omega) \end{align}

The difference here is that $$u_h$$ is generally not continuous differentiable on the whole $$\Omega$$, but only on the elements. Therefore you are only allowed to apply the div-theorem element wise. Note that in general:

\begin{align} \sum_{T\in \mathcal{T}} \int_{\partial T}\nabla u_h \cdot n v_h \ ds \neq \int_{\Gamma} \nabla u_h\cdot n v_h \ ds = 0 \end{align}

So the normal-flux of the solution $$u_h$$ dosent need to cancel on the faces inside $$\Omega$$. Nevertheless $$u_h$$ will give you an approximation of the original solution $$u$$.

On the other side for a FVM you have by applying the div-theorem the weak form:

\begin{align} -\sum_{T\in \mathcal{T}}\int_{\partial T} \nabla u \cdot n \ ds= -\int_\Gamma \nabla u \cdot n \ ds = \sum_{T\in \mathcal{T}}\int_T f \ dx \end{align}

And further even in the discrete setting the contributions inside the domain cancel since

\begin{align} -\sum_{T\in \mathcal{T}}\int_{\partial T} \nabla u_h \cdot n \ ds = \sum_{T\in \mathcal{T}}\int_T f \ dx = -\int_\Gamma \nabla u_h \cdot n \ ds \end{align}

This is why we call FVM often conservative by construction.

• How come $u$ gets to be continuously differentiable but $u_h$ doesn't if $V_h$ is a subset of $V=H^1$ ? Oct 17, 2021 at 3:39
• I say this because if you assume that $f\in L_2$ then $\Delta u\in L_2$. Moreover if you have a convex domain then by the shift theorem $u\in H^2$. So this is why I say here that $u$ is continuously differentiable. Now if you take $V_h=\{v_h\in H^1(\Omega):v_h|T=\mathcal{P^k}(T)\ \forall T\}$ you see that only the first derivative exists in a weak sense and can jump between elements.
– Pepe
Oct 17, 2021 at 10:11

The face integrals don't show up when using continuous basis functions, and the weak form that you wrote down indeed assumes that the basis functions are continuous. The finite volume method is a special case of the finite element method when you allow discontinuous basis functions instead of only continuous ones. For a discontinuous Galerkin (or DG) discretization of the Poisson equation there are indeed face terms. I wrote a bit about the derivation and motivation for this here if you want to read more. But the short answer is that the weak form of the symmetric interior penalty DG method is that

$$\sum_K\int_K\left(\nabla u_h\cdot\nabla v_h - f\cdot v_h\right)dx + \sum_E\int_E\left([u_h]\left[\frac{\partial v_h}{\partial n}\right] + \left[\frac{\partial u_h}{\partial n}\right][v_h] + \frac{\gamma}{h}[u_h][v_h]\right)dS = 0$$

for all test functions $$v_h$$, where $$K$$ are all the cells of the mesh, $$E$$ the faces, $$\gamma$$ is a penalty parameter that can be calculated explicitly in terms of the element degree and mesh quality, and $$[\cdot]$$ denotes the jump of a quantity across cell boundaries.

• It is still weird for me that the weak form doesn't have something like that built in. Is it because we are no long in the $H^1$ space? If that is the case how can we tell that the $H^1$ space is the space that the solution MUST be in. I guess so. I know it is the natural space but if you have any insights I'd like to know.. But isn't the space of local support functions discontinuous? Oct 13, 2021 at 1:27
• You're right, it is weird! The reason is that DG discretizations are non-conforming -- the basis functions don't live in $H^1$, you're approximating a function in that space "from outside" in a manner of speaking. Oct 13, 2021 at 3:59
• Let me add that the solution space depends on the weak formulation. In your weak formulation there is an integral $\int_\Omega \nabla u^2 \, \mathrm{d}x$ (after choosing $v=u$) which is finite if and only if $u \in H^1(\Omega)$. So that's why the solution MUST be in $H^1$. However, in DG your discrete solution is not (necessarily, most of the time) in $H^1$ but in $L^2$. This discrete function can still approximate a function in $H^1$ though.
– knl
Oct 14, 2021 at 7:26