1
$\begingroup$

I would like to ask how is the evaluation of integrals over inter-element interfaces implemented in a typical DG code. I can think of two basic approaches (I assume 2D mesh here):

1) Perform 1D interpolation of the interface data seen from 'left' and 'right' of each edge and numerically integrate. For this, I don't need to evaluate the element shape functions in the elements adjacent to the edge. This will probably work fine for Lagrange shape functions (element interpolant restricted to edge and 1D edge interpolant should be identical).

2) Interpret the 1D edge quadrature rule as a quadrature with points placed on the boundary of each element and use the element shape functions to compute edge integrals.

Option 1) does not require extra storage for jump/average terms: evaluation on each edge is performed on the fly and I can directly add boundary integrals to the weak form. Option 2) will require a lot of extra memory allocation for interface integrals: jump and average terms are only complete after I visit both left and right neighbour of each edge and temporary data evaluated at quadrature points have to wait somewhere before the interface integrals are assembled. I could also loop over all edges and obtain directly the left-and right interpolant from the element shape functions, but that would mean that I re-evaluate the element shape functions every time I visit one of its edges (computational overhead).

I think that the disadvantage of 1) is that it won't be possible to evaluate shape function gradients on element interfaces for higher-order equations. In addition, it would probably not work with some modal bases (all modes from the element interior are needed in the expansion).

I cannot decide which approach is better (or if they are good at all). Maybe there's a more clever way I did not think of. Can you share your opinion?

Improve this question
$\endgroup$
1
  • 1
    $\begingroup$ Are you using a nodal or modal DG method? With nodal the choice is pretty clear: you already have the values for interfacial fluxes at the quadrature points. I'm not fully understanding your question -- for DG, the surface integrals in the weak form always depend on both sides of the face. As for the higher order equations, that's an open research question. See the Bassi & Rebay paper (BR2) for one common method for 2nd derivatives. $\endgroup$
    – Aurelius
    Commented Jun 25, 2014 at 19:09

1 Answer 1

Reset to default
1
$\begingroup$

In deal.II (see http://www.dealii.org/), we use your option 2. The reason is that you need to evaluate not just the shape functions that are defined along the interface, but you also need to evaluate all of the other shape functions defined on that cell (even if their interpolation point -- in the case of Lagrange elements -- is not on the given interface: while they will be zero along the interface, their gradient is not).

You can see how this is implemented in the step-12 tutorial program of deal.II: http://www.dealii.org/developer/doxygen/deal.II/step_12.html

Improve this answer
$\endgroup$
2
  • $\begingroup$ Ok, thank you for your answer and for the link. I still don't understand where (in the element) would you evaluate the shape function gradients (if you needed them) in integrate_face_term: on the common face of cell1 and cell2 only? This means that when you process the interface of cell1 and say cell3, you will again have to evaluate all shape function gradients in cell1, just on different face. Or do you directly compute the gradients on all faces of cell1 and store temporary data somewhere until you visit cell3? Wouldn't this require a lot of memory? $\endgroup$
    – Martin Vymazal
    Commented Jun 26, 2014 at 10:40
  • $\begingroup$ This happens in the FEFaceValues objects: They evaluate the shape functions (and other things, as necessary) at quadrature points on the face you select. Once you move on to another face, the location of the quadrature points changes and so the evaluation has to be repeated. $\endgroup$
    – Wolfgang Bangerth
    Commented Jun 29, 2014 at 12:30

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.