Gradient-jump penalty term in FEM

Question

I am slightly confused regarding the meaning of the $i-th$ gradient-jump term $[\nabla \phi_i]$ in the context of finite element methods, used in the assembly of the stiffness matrix (an example with deal.ii). On face $F$ with unit normal $\nu$, the jump term is defined as $$ [\nabla \phi ] = \nu \cdot (\nabla\phi^+-\nabla\phi^-) $$ For example, consider the case of a Lagrangian, bilinear quadrilateral element, where $i \in \{0,1,2,3\}$. We want to compute the gradient-jump penalty term on face $F$ \begin{gather} \int_F [\nabla u][\nabla \phi] \tag{$\star$} \end{gather} But we see that, since $u = \sum_i \phi_i u_i$ and $u_0^+ \equiv u_3^-$, $u_1^+ \equiv u_2^-$ (based on the numbering depicted in the figure below) \begin{align} \nabla u^+ - \nabla u^- &= (\nabla \phi_0^+-\nabla \phi_3^-)\ u_0^+ + (\nabla \phi_1^+ - \nabla \phi_2^-)\ u^+_1 \\ &+ \nabla \phi_2^+\ u_2^+ - \nabla \phi_1^- u_1^- + \nabla \phi_3^+\ u_3^+ - \nabla \phi_0^-\ u_0^- \end{align}

So I observe that for the shared vertices the terms can be grouped and the calculations can be made on the $+$ element but for the non-shared vertices there are terms for both elements ($+,-$).

That being said, what is the meaning of $[\nabla \phi_i]$ if $(\star)$ cannot be written in the form $$ \sum_j \int_F [\nabla \phi_i] [\nabla \phi_j]\ u_j \quad ? $$

Note.

I can see that in deal.ii there is a class named FEInterfaceValues and probably what I am asking is what its function dofmap is actually doing, but I don't understand the implementation.

Answer 1

You go the grouping in the term $\nabla u^+-\nabla u^-$ mostly right but not quite. To see this, remember that the functions $\phi_i^{\pm}$ are discontinuous, and that consequently $u_0^+$ is not necessarily equal to $u_3^-$ -- if they were, you'd end up with a function that is continuous at that point. Instead, you need one degree of freedom per adjacent cell at each vertex.

The better way to write things is always to say that the shape functions form a basis. If you want to stick with the mesh with the two cells, then you will have eight shape functions $\phi_i^\pm$ (for $i=0\ldots 3$) and you need eight coefficients $u_i^\pm$. Your solution $u_h$ is then simply given by $u_h=\sum_{i=0}^3 \phi_i^- u_i^- + \sum_{i=0}^3 \phi_i^+ u_i^+$.

From this you can then see what the jump in the gradient is going to be: $$ [\nabla u] = \left[\nabla\left(\sum_{i=0}^3 \phi_i^- u_i^- + \sum_{i=0}^3 \phi_i^+ u_i^+\right)\right] \\ = \sum_{i=0}^3 [ \nabla \phi_i^-] u_i^- + \sum_{i=0}^3 [\nabla\phi_i^+] u_i^+. $$ So your expansion in the last formula of the question is pretty much exactly right if you re-interpret the $u_j$ as my $u_i^\pm$.

What FEInterfaceValues does is pretty much the following: For a face, it gathers the indices of all of the shape functions that live on the neighboring cells. In your case, that's 8 shape functions, and so when assembling these jump terms along the face, you'll have to compute an $8\times 8$ matrix of jump contributions to the global matrix. But if you had, for example, continuous elements, there would only be 6 (unique) degrees of freedom adjacent, and so the dofmap variable would only contain 6 entries (of the global indices of the corresponding degrees of freedom).