# Projection of user expression onto DG1

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(),

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()
``````

Dolfin 2019.2, Ubuntu LTS 20.04, PPA

The `"ufl"` implementation evaluates `ff` at quadrature points on the interiors of elements during the `project`ion, 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.

What is the reason for such behaviour? Can we force Fenics evaluates â€śnumericâ€ť at quadrature points?

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:

# Function space with DoFs at quadrature points:

QE = df.FiniteElement(family=fe_space, cell=mesh.ufl_cell(),

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

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

Hi David,
thank you. It works!