Projection of user expression onto DG1

1

Hi,
I want to approximate a discontinuous function in DG1 space by L2 projection. Why does the choice of the user or ufl expression cause a difference in the approximation? Please, see the following MWE:

import dolfin as df
import matplotlib.pylab as plt
import numpy as np

class ScalarFunction(df.UserExpression):
    def __init__(self, fun_to_eval, **kwargs):
        super().__init__(**kwargs)
        self.fun_to_eval = fun_to_eval

    def eval(self, value, x):
        value[0] = self.fun_to_eval(x)
    
    def value_shape(self):
        return ()


def fun_to_eval(x):
    
    if x >= 0.3:
        return x**2
    else:
        return x-1.
    
if __name__ == '__main__':
    mesh = df.IntervalMesh(10, 0.0, 1.0)
    x = df.MeshCoordinates(mesh)

    fe_space = 'DG'
    fe_degree = 1
    expr = 'numeric'
    
    QE = df.FiniteElement(family=fe_space, cell=mesh.ufl_cell(),
                          degree=fe_degree, quad_scheme="default")

    V = df.FunctionSpace(mesh, QE)
    
    if expr=='ufl':
        # ufl expression
        ff = df.conditional(df.ge(x[0], 0.3), x[0]**2, x[0]-1.)
    elif expr=='numeric':
        # numerical expression
        ff = ScalarFunction(fun_to_eval, element=QE)

    f = df.project(ff, V)
    
    xx = np.linspace(0, 1, 50)
    fx = np.zeros(len(xx))
    fx_p = np.zeros(len(xx))
    for i in range(len(xx)):
        fx[i] = f(xx[i])
        fx_p[i] = fun_to_eval(xx[i])
        
        
    plt.plot(xx, fx_p, label='analytical')
    plt.plot(xx, fx, 'o', label='FE')
    plt.title(expr + '-' + fe_space + str(fe_degree))
    plt.legend()
    plt.show()

results
results3217Ă—1100 161 KB

Dolfin 2019.2, Ubuntu LTS 20.04, PPA

2

The "ufl" implementation evaluates ff at quadrature points on the interiors of elements during the projection, so the values are different on each side of the discontinuity at x=0.3. The "numeric" implementation is first interpolating ScalarFunction at nodes of the space V before projecting. The boundary nodes of the DG elements to either side of the discontinuity will have the same location, and will therefore interpolate to the same value, causing ff to be continuous.

3

What is the reason for such behaviour? Can we force Fenics evaluates “numeric” at quadrature points?

4

You can pass an element of type "Quadrature" instead of "DG" to the constructor for ScalarFunction to interpolate the UserExpression at quadrature points:

import dolfin as df
import matplotlib.pylab as plt
import numpy as np

class ScalarFunction(df.UserExpression):
    def __init__(self, fun_to_eval, **kwargs):
        super().__init__(**kwargs)
        self.fun_to_eval = fun_to_eval

    def eval(self, value, x):
        value[0] = self.fun_to_eval(x)
    
    def value_shape(self):
        return ()

# Re-define `project` function to allow
# explicit specification of an integration
# measure (to ensure consistent quadrature degree
# between quadrature space and projection assembly).
def custom_project(f,V,dx):
    Pf = df.TrialFunction(V)
    v = df.TestFunction(V)
    Pf_h = df.Function(V)
    df.solve(df.inner(Pf,v)*dx==df.inner(f,v)*dx,Pf_h)
    return Pf_h

def fun_to_eval(x):
    
    if x >= 0.3:
        return x**2
    else:
        return x-1.
    
if __name__ == '__main__':
    mesh = df.IntervalMesh(10, 0.0, 1.0)
    x = df.MeshCoordinates(mesh)

    fe_space = 'DG'
    fe_degree = 1
    expr = 'numeric'

    # Must use the same quadrature degree for everything, so
    # quantities in quadrature spaces are well-defined:
    QUAD_DEG = 2
    dx_deg = df.dx(metadata={"quadrature_degree":QUAD_DEG})
    
    # Function space with DoFs at quadrature points:
    QE_quad = df.FiniteElement("Quadrature", mesh.ufl_cell(),
                               degree=QUAD_DEG,
                               quad_scheme="default")
    
    QE = df.FiniteElement(family=fe_space, cell=mesh.ufl_cell(),
                          degree=fe_degree, quad_scheme="default")

    V = df.FunctionSpace(mesh, QE)
    V_quad = df.FunctionSpace(mesh, QE_quad)
    
    if expr=='ufl':
        # ufl expression
        ff = df.conditional(df.ge(x[0], 0.3), x[0]**2, x[0]-1.)
    elif expr=='numeric':
        # numerical expression
        
        # Interpolate on quadrature element instead of DG:
        ff = ScalarFunction(fun_to_eval, element=QE_quad)

    f = custom_project(ff, V, dx_deg)
    
    xx = np.linspace(0, 1, 50)
    fx = np.zeros(len(xx))
    fx_p = np.zeros(len(xx))
    for i in range(len(xx)):
        fx[i] = f(xx[i])
        fx_p[i] = fun_to_eval(xx[i])
        
    plt.plot(xx, fx_p, label='analytical')
    plt.plot(xx, fx, 'o', label='FE')
    plt.title(expr + '-' + fe_space + str(fe_degree))
    plt.legend()
    plt.show()

"Quadrature" function spaces have degrees of freedom at integration points of a specified quadrature rule. Some notes for use of these spaces:

  • If a Function of type "Quadrature" appears in a Form, that Form must be integrated with a Measure of the same quadrature scheme/degree. As in the above example, this typically requires you to manually override FEniCS’s default behavior of automatically estimating an appropriate quadrature degree.
  • If a Function of type "Quadrature" is the unknown in a variational problem (which is not the case here, since ff is still projected into the DG space V, not V_quad), you need to circumvent a bug in the default form compiler (UFLACS) by specifying the (deprecated) "quadrature" form representation, as demonstrated in @bleyerj’s tutorial here (where he uses "Quadrature" elements to represent local history variables in plasticity).
4 Likes
5

Hi David,
thank you. It works!