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 ? $$


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.

Schematic of the FEM stencil

  • 3
    $\begingroup$ I think I'm reasonably well positioned to answer the question :-) But I don't follow how you arrive at $\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^- $. Can you explain this formula? $\endgroup$ – Wolfgang Bangerth Mar 3 at 13:54
  • $\begingroup$ Separately, are you assuming that your finite element is continuous or discontinuous? $\endgroup$ – Wolfgang Bangerth Mar 3 at 13:55
  • $\begingroup$ @WolfgangBangerth Thank you for your reply. I updated the post to explain the formula. I am assuming bilinear quadrilaterals, but added the DG flag since similar terms appear in DG-FEM. $\endgroup$ – ares Mar 3 at 15:52
  • $\begingroup$ So does $\phi_3^-$ live on the bottom cell, and at the same location as $\phi_0^+$? $\endgroup$ – Wolfgang Bangerth Mar 3 at 20:04
  • $\begingroup$ @WolfgangBangerth Exactly. Sorry if that was not clear. $\endgroup$ – ares Mar 3 at 22:35

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).

| cite | improve this answer | |
  • $\begingroup$ Thank you for the detailed answer. In my question I assumed continuous shape functions, hence the grouping. Could you please define the terms $[\nabla \phi_i^+]$ and $[\nabla \phi_i^-]$? To me it looks like $[\nabla \phi_i^+] \equiv \nu \cdot \nabla \phi_i^+$, since $\nabla \phi_i^+ \equiv 0$ if you approach from the `$-$' FE. $\endgroup$ – ares Mar 4 at 12:32
  • $\begingroup$ Thank you for the clarification regarding the deal.ii functions. Currently I am writing my own code but I may have to migrate to deal.ii for bigger problems, so it's good to know it's already implemented. $\endgroup$ – ares Mar 4 at 12:38
  • $\begingroup$ Yes, your expansion of what the jump of a shape function is is correct in the case of a discontinuous shape function. $\endgroup$ – Wolfgang Bangerth Mar 5 at 0:00
  • $\begingroup$ Shouldn't $[\nabla \phi_i^+] = \nu \cdot \nabla \phi_i^+$ and $[\nabla \phi_i^-] = \nu \cdot \nabla \phi_i^-$ hold for any shape function since $\phi^+, \ \phi^-$ have support only on $+$ and $-$ FE, respectively? Why would it have to be a discontinuous one? $\endgroup$ – ares Mar 5 at 0:15
  • 1
    $\begingroup$ No, for continuous functions, $\phi_0^+$ and $\phi_3^-$ would denote the same shape function. Similarly, $u_0^+=u_3^-$ are the same unknowns. $\endgroup$ – Wolfgang Bangerth Mar 5 at 1:37

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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