I am currently learning nodal DG methods, primarily through the book by Warburton, and am a bit confused on how to handle surface integrals using straight edged elements. On page 187 (and on page 214) of Warburton's book, we have

$$ \int_{\partial E_{k}} \mathbf{n} \cdot f_{h} \ell_{i}^{k}(\mathbf{x}) d\mathbf{x}. $$

Recall that $f_{h} = \sum_{j=1}^{N} f_{j}(x,t)\ell_{j}(\mathbf{x})$ and $\ell_{i}(\mathbf{x})$ is our basis function, and so,

$$ \int_{\partial E_{k}} \mathbf{n} \cdot f_{h} \ell_{i}^{k}(\mathbf{x}) d\mathbf{x} = \int_{\partial E_{k}} \mathbf{n} \cdot \left[ \sum_{j=1}^{N} f_{j}(x,t)\ell_{j}(\mathbf{x}) \right]\ell_{i}^{k}(\mathbf{x}) d\mathbf{x}\\ = \sum_{j=1}^{N} \mathbf{n} \cdot f_{j}(x,t) \int_{\partial E_{k}} \ell{j}(\mathbf{x}) \ell_{i}(\mathbf{x}) d \mathbf{x}\\ = \sum_{j=1}^{N} \mathbf{n} \cdot f_{j}(x,t) M^{k}_{ij}. $$ Now, to discretize the surface mass matrix, we must use a quadrature rule on each face (since I am working in 2D it's a line integral along each edge). With that said, since we have an orthonormal basis, we can compute the mass matrix analytically. For volume integrals, Warburton shows that $M^{k}_{ij} = (\mathcal{V} \mathcal{V}^{T})^{-1} J,$ where $J$ is the metric Jacobian. So intuitively, the same would be for the mass matrix, however, the Vandermonde matrix, $\mathcal{V}$ would be constructed so its on the local face points $(r,s)$, where $(r,s) \in [-1,1]$. Now, looking at their code (https://github.com/tcew/nodal-dg/blob/master/Codes1.1/CFD2D/CurvedEulerRHS2D.m), we see that the surface integral term has no surface mass matrix. There is an inverse mass matrix applied to the numerical flux but it is the one from the time-derivative. I was curious if someone could point out why there is no surface mass matrix in their formulation.

  • $\begingroup$ I assume that $\ell_j$ is a basis function, but it won't hurt if you clarify your notation. $\endgroup$
    – nicoguaro
    Commented Jun 5, 2021 at 14:25
  • $\begingroup$ For what I understand you I think that you want to interpolate your flux using the same basis functions, that's why you get a mass matrix. But, that's not the only way to solve that integral. You could, for example, use numerical integration and evaluate the function in a discrete set of points. Also, you could assume that the function is constant in your element and use the centroid value. $\endgroup$
    – nicoguaro
    Commented Jun 5, 2021 at 14:28
  • $\begingroup$ I edited my post to specify what l_{j}(x) is. I understand you want to numerically integrate the mass matrix at the quadrature points, however, what I don't understand is why it seems to be 1 in this case (specifically line 56 of the code). In their example code, they interpolate the f_{j} and n to the edges using Vf (Vf = Vg*inv(V), where Vg is the Vandermonde matrix evaluated at the quadrature points). However, why would that remove the need to directly evaluate the mass matrix as the edges? Sorry if I am completely missing something. Wouldnt you interpolate f_{j}, n, and M to the edges? $\endgroup$
    – David1998
    Commented Jun 5, 2021 at 14:48
  • $\begingroup$ You don't need a mass matrix for surface terms. That appears only if you interpolate $f$. I haven't taken a look at the code. $\endgroup$
    – nicoguaro
    Commented Jun 5, 2021 at 15:25

1 Answer 1


In 2D your volume part consists of a double integral, whereas the surface part is a standard line integral. Generally this is described with tensorial notations.

If you consider the Cartesian case, the volume part results in a mass matrix similar to

$M_{\Omega}\equiv \mathbf{M}_{\text{1}} \otimes \mathbf{M}_{\text{1}}$,

whereas the surface part results in a mass matrix similar to

$M_{\partial \Omega}^{\xi} \equiv \mathbf{M}_{\text{1}} \otimes \mathbf{I}_1\\ M_{\partial \Omega}^{\eta} \equiv \mathbf{I}_1 \otimes \mathbf{M}_{\text{1}}$.

Here $\mathbf{I}_1$ is the identity matrix.

If you invert the volume part and apply it on $\xi$-faces you will get something like

$=\left( \mathbf{M}_{\text{1}} \otimes \mathbf{M}_{\text{1}}\right)^{-1}(\mathbf{M}_{\text{1}} \otimes \mathbf{I}_1) \\ = (\mathbf{M}_{\text{1}} ^{-1} \otimes \mathbf{M}_{\text{1}} ^{-1}) (\mathbf{M}_{\text{1}} \otimes \mathbf{I}_1) \\ = (\mathbf{M}_{\text{1}} ^{-1} \mathbf{M}_{\text{1}} ) \otimes (\mathbf{M}_{\text{1}}^{-1} \mathbf{I}_1)\\ = \mathbf{I}_{\text{1}} \otimes \mathbf{M}_{\text{1}}^{-1}$

or similar on $\eta$-faces

$=\left( \mathbf{M}_{\text{1}} \otimes \mathbf{M}_{\text{1}}\right)^{-1}(\mathbf{I}_{\text{1}} \otimes \mathbf{M}_1) \\ = (\mathbf{M}_{\text{1}} ^{-1} \otimes \mathbf{M}_{\text{1}} ^{-1}) (\mathbf{I}_{\text{1}} \otimes \mathbf{M}_1) \\ = (\mathbf{M}_{\text{1}} ^{-1} \mathbf{I}_{\text{1}} ) \otimes (\mathbf{M}_{\text{1}}^{-1} \mathbf{M}_1)\\ = \mathbf{M}_{\text{1}}^{-1} \otimes \mathbf{I}_{\text{1}}.$

Summerizing: I think you are missing the fact, that both mass matrices (volume and surface) cancel, resulting in a 1D inverse mass matrix. Note that numerical fluxes on $\eta$-faces are projected into the volume in $\xi$-direction, whereas fluxes on $\xi$-faces are projected into the volume in $\eta$-direction.



Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.