Pointwise expressions

[garth-wells]

From @blechta

From time to time there appears desire for pointwise expressions, i.e., expressions which get evaluated in forms at quadrature points rather than interpolated to local FE space (typically Lagrange). Absurdly, FEniCS has support for pointwise expressions for a long time through quadrature element, see https://bitbucket.org/snippets/blechta/Rezng6.

I argue that implementation through quadrature element is a right approach. It just fits the pipeline:

• UFL: representation, preprocessing, including degree estimation,
• tabulate_tensor-interface: no need to have special kernel interface for pointwise expression; it is passed as expansion coefficients of quadrature element; remind that we want simple C replacement of UFC so this is a right approach,
• C++ DOLFIN: nothing new is needed on C++ side; dolfin::Expression
• Python DOLFIN: it is possible to have pointwise expression which is cpp.Expression and ufl.Coefficient on quadrature element but it is arguable how this multiple inheritance is confusing and sustainable, see Separate cpp.Expression and ufl.Coefficient #1.

Taking the assumption that the quadrature element approach is a right approach, the following text is relevant. Range of approaches is available, providing more or less magic functionality and difficulty of implementation. The dilemma comes from the fact that

e = PointwiseExpression("cos(x[0])")
e*u*v*dx

does not have any sense as quadrature degree is not specified. It can be done in many several ways and all of them have advantages and disadvantages.

1. Soft, trivial approach: provide documentation and coverage for existing functionality. With that users and developers would be aware how to specify pointwise expressions:

qe = FiniteElement('Quadrature', mesh.ufl_cell(), 6, quad_scheme="default")
my_pointwise_expression = MyExpr("x[0]*x[0]", element=qe)

# Compatible quadrature specification is needed
integral = my_pointwise_expression*other_coefficients*dx(metadata={'quadrature_degree': 6})

The advantage is that there is no magic functionality (no need for some tricky UFL preprocessing as element is fully specified from the beginning). Disadvantage is that some people might find it pretty ugly and discouraging from usage.

2. Less magic user friendly factory. It turns out it is easy to support:

my_pointwise_expression = create_pointwise_expression("x[0]*x[0]", quadrature_degree=6)
integral = my_pointwise_expression*other_coefficients*dx

Here user calls the factory function with compiled or Python expression and total quadrature degree in the integral where they intends to use the expression. Arguably this is already too much from mathematical perspective because it hides quadrature details from the measure into the expression. (Less invasive approach of still requiring degree specification in measure is possible.)

3. Very magic, super-friendly factory. It turns out that it is possible to implement:

# Expression contributes to total quadrature by 2
my_pointwise_expression = create_pointwise_expression("x[0]*x[0]", degree=2)
integral = my_pointwise_expression*other_coefficients*dx  # quad degree is 2 + estimated degree of other_coefficients

# Expression does  not have degree; quadrature degree specified in measure
my_pointwise_expression = create_pointwise_expression("x[0]*x[0]")
integral = my_pointwise_expression*other_coefficients*dx(metadata={'quadrature_degree': 6})

Both of these can be supported simultaneously. There is a working sketch of this. Unfortunately it turns out that this is very problematic if the form consists of more integrals (possibly having different quadrature degree). To support this properly, major change in UFL form preprocessing and form compilers would be needed (coefficient replace maps would need to be per integral rather than per form).

Approach 2. is easily possible now, but might require the same UFL preprocessing overhaul if mixed cells support is wanted.