Projection of user expression onto DG1

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).
3 Likes