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

variational formulation
1

Hi, I am trying to implement a fictitious domain method for a fluid simulation. This requires an evaluation of a penalty function which is discontinuous within the element. In theory, the penalty function should be evaluated on each Gauss quadrature points during assembly.

It seems like FEniCS evaluate UserExpression at the nodes of the element and use interpolation to get the value of the the penalty function at each Gauss quadrature points when assemble is called. Is there a way to make FEniCS evaluate the UserExpression on each Gauss quadrature point during assembly or if there are any workarounds for this?

Attached here is a small toy problem created to test the evaluation process. The UserExpression defines the area in which the penalty function is applied. Here it returns 1 if the point is inside the circle and 0 of the point is outside of the circle. The script integrates (1-circle)*dx which on paper should give (1-pi r^2) where r is the circle radius. If FEniCS evaluate the UserExpression at each Gauss quadrature points, then as quadrature degree increases, the integration should approach (1-pi r^2) which is what I’d like to achieve in my simulation.

import fenics as fe
import matplotlib.pyplot as plt
import numpy as np

# -- quadrature degree, expression degree, and function space degree
quadDeg = 1
exprDeg = 2
funSpaceDeg = 1
fe.parameters["form_compiler"]["quadrature_degree"] = quadDeg

# -- Inside-Outside userexpression
class InsideCircle(fe.UserExpression):
    def __init__(self,posx,posy,r,**kwargs):
        self.m_posx = posx
        self.m_posy = posy
        self.m_r = r
        super().__init__(**kwargs)
    def eval(self, value, x):
        # -- Return 1 if x is inside the and 0 if x is outside of the circle
        if ( (x[0]-self.m_posx)**2 + (x[1]-self.m_posy)**2 ) < self.m_r**2:
            value[0] = 1
        else:
            value[0] = 0
        print( 'Evaluating at position: (%f,%f). Value: %d' % (x[0],x[1], value[0]) )
        return value
    def value_shape(self):
        return ()

# -- Define circle position and redius
posx = 1/2
posy = 1/2
r = 0.2
circle = InsideCircle(posx,posy,r,degree=exprDeg)

# -- Integrate 
mesh = fe.UnitSquareMesh(1,1)
V = fe.FunctionSpace(mesh,'Lagrange',funSpaceDeg)
one = fe.Expression('1',degree=funSpaceDeg)
one = fe.interpolate(one,V)
area=fe.assemble((one-circle)*fe.dx)
print('Computed area: %f' % (area))
print('Actual area: %f' % (1-np.pi*r**2))

UnitMesh3

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