General concepts » Interfaces, traces, and discontinuous Galerkin

Working on facets: traces, jumps and averages, interface tagging, and DG formulations.

Introduction

Many formulations live on facets rather than cells: Neumann data on the boundary, transmission conditions on material interfaces, and the entire discontinuous Galerkin (DG) family. This guide collects the facet-side vocabulary of the form language and its correct use.

The facet integral family

ClassRanges overTypical use
BoundaryIntegralExterior (codimension-1 boundary) facetsNeumann/Robin terms
InterfaceIntegralInterior facetsDG coupling, transmission terms
FaceIntegralAll facetsTerms defined uniformly on faces

All accept attribute restriction via .over(attr, ...), so a term can target a tagged subset of facets. Normals are provided by BoundaryNormal (outward on the boundary) and FaceNormal (a fixed orientation convention on interior faces — sign-sensitive terms must use it consistently on both sides).

Traces: which side is meant?

An $ H^1 $ function has a well-defined boundary trace, but on an interior* facet a discontinuous quantity (a P0 function, a gradient of a P1 function, any DG field) has two one-sided values. Rodin requires you to disambiguate:

  • f.traceOf(attr) selects the trace seen from the cells with attribute attr — the tool for transmission problems where a coefficient jumps across a material interface;
  • if evaluation on an interior facet is attempted without a determined side, the library raises UndeterminedTraceDomainException — this exception almost always means a missing traceOf on a coefficient in an interface term.

Interfaces themselves are just tagged facets: use mesh.trace({{attr1, attr2}, ifaceAttr}) to label the facets between two cell attributes, after computing the facet-cell connectivity (see Connectivity).

Jump and average

The two facet operators of DG methods, defined for a quantity $ w $ with one-sided values $ w^+, w^- $ :

\[ [\![ w ]\!] = w^+ - w^-, \qquad \{ w \} = \tfrac12 (w^+ + w^-), \]

available as Jump and Average:

auto jf = Jump(f);       // f: any function with facet traces
auto af = Average(f);

Sanity properties worth internalizing (they are what the unit tests assert): the jump of any continuous function — a constant, or a conforming P1/H1 grid function — is zero on every interior facet; the average of a continuous function equals its trace.

Interior penalty DG in the form language

The symmetric interior penalty (SIPG) discretization of $ -\Delta u = f $ reads

\[ \sum_K \int_K \nabla u \cdot \nabla v \;-\; \int_{\mathcal F_i} \{\nabla u\}\cdot n \,[\![ v ]\!] \;-\; \int_{\mathcal F_i} \{\nabla v\}\cdot n \,[\![ u ]\!] \;+\; \int_{\mathcal F_i} \frac{\eta}{h} [\![ u ]\!]\,[\![ v ]\!] \;=\; \int_\Omega f v , \]

written with Integral for the cell term and InterfaceIntegral for the three facet terms (plus BoundaryIntegral analogues to impose Dirichlet data weakly, Nitsche-style). Structural facts the mathematics imposes:

  • the penalty scaling $ \eta / h $ with $ \eta $ large enough (growing with the polynomial degree) is what makes the form coercive — omitting or under-sizing it is the standard DG instability, and it manifests as an indefinite system, not a crash;
  • the two consistency terms must be transposes of each other for symmetry (SIPG); flipping one sign gives the nonsymmetric NIPG variant with different penalty requirements;
  • assembled DG systems are not SPD in general — see Solvers for appropriate solvers.

See examples/DG for a starting point, and the Jump/Average unit tests (tests/unit/Rodin/Variational/JumpTest.cpp, AverageTest.cpp) for the precise operational semantics.

See also

« Solvers | General concepts | Shape and topology optimization »