Skip to content

Precise semantics of `LocalFiniteInterpolation::interpolate()`

(This issue is about the local interpolation interface and should thus probably be discussed in dune-localfunctions. However, the question mainly comes up when using this in a global context. Hence I consider it more practical to discuss this here first. The result may serve as a well founded proposal with an example implementation in dune-functions.)

For a LocalFiniteElement (LFE) its LocalInterpolation::interpolate(f, coeff) should evaluate the dual basis functions in the sense of Ciarlet at f and write the computed coefficients. Since we're doing everything locally, f has to live on the reference element. Assuming that a globally defined function g should be interpolated, the local function f satisfies f(x)=g(Phi(x)) where Phi is the element transformation. When using Hermite elements, interpolate() has to evaluate derivatives of f the question shows up, which type of derivatives of f should be used. There's two possibilities:

  • (A) interpolate() should use Df(x) = Dg(Phi(x))*DPhi(x).
  • (B) interpolate() should use Dg(Phi(x))

Arguments for (A):

  • Conceptually, the local finite element does not know about some global transformation and should expect everything to be local and hence simply use the derivative of f.
  • (A) ensures the local bi-orthogonality relation between primal and dual basis functions in the sense of Ciarlet on the reference element.

Arguments for (B):

  • (B) is what we decided on for our grid-function interface: localFunction(g) is not the transformation of f onto some element but encodes its restriction onto the element which happens to be defined in terms of local coordinates.
  • (B) allows to pass a g_local = localFunction(g) directly to interpolate.
  • (A) would require to build a wrapper f around g_local that where f(x)=g_local(x) but derivative(f) = derivative(g_local(x))*DPhi(x).
  • For a Hermite-element the local finite element does know about the transformation, because its local basis functions depend on it.
  • When implementing the interpolate() it is often very natural to implement the evaluation of (B). Using Df as in (A) may require to internally transform this back to (B).
  • The internal transformation means that we do a lot of useless work for (A): First we build (outside) the wrapper f that multiplies with DPhi(x). This is passed to interpolate(f, ...) that internally multiplies with DPhi(x)^{-1}.
  • Even worse, when passing a DiscreteGlobalBasisFunction to the global interpolate() method, then derivative(localFunction(g)) has to transform the local basis Jacobians with DPhi(x)^{-1} which leads to three transformations in a row, where two cancel out.

That's a lot of arguments for (B), but still I consider the consistency argument for (A) to be very valid. Maybe we a real solution would be:

  • (C) Introduce a global finite element interface (cf. in dune-localfunctions). This follow our localFunction decision and represent the actual finite element on some real geometry, but happen to be implemented in terms of global coordinates. Hence this should naturally assume (B) and also transform the derivatives in evaluateJacobian(), such that the bi-orthogonality is preserved.

(C) has the disadvantage, that it requires to do significant additions to our interfaces, especially we would have two concurrent interfaces in our leaf nodes a local and a global FE. Then we would have to decide which of them are optional, and how our algorithms should decide which one to use. In contrast to this (A) and (B) are rather light-weight, because they have no influence an non-Hermite elements. Hence I would opt to first do (A) or (B) to get started.

BTW: !397 nicely reflects the fact that (B) is very natural in terms of the implementation.

Edited by Carsten Gräser