How to evaluate a UserExpression on each gauss-quadrature point during assembly?

variational formulation
2

It is possible to pass an element type instead of a degree to the UserExpression's constructor. One of the element types in FEniCS is “Quadrature”, which has degrees of freedom at each quadrature point. You can use this in your code as follows:

#circle = InsideCircle(posx,posy,r,degree=exprDeg)

# -- Integrate
mesh = fe.UnitSquareMesh(1,1)

QE = fe.FiniteElement(family="Quadrature", cell=mesh.ufl_cell(),
                      degree=quadDeg, quad_scheme="default")
circle = InsideCircle(posx,posy,r,element=QE)
fe.dx = fe.dx(metadata={"quadrature_degree":quadDeg})

This appears to converge to the correct answer as quadDeg is increased, but I’ll add a few notes:

  1. Quadrature elements are somewhat buggy in FEniCS. Using quadrature function spaces for unknown functions in variational problems is still not well-supported (although there are some workarounds).

  2. Gaussian quadrature assumes that the integrand is well-approximated by a polynomial interpolant. Interpolating a discontinuous function (such as the characteristic function of a circle) with a high-order polynomial gives a highly-oscillatory result.

1 Like