Hi,
I want to use the custom quadrature rule patch proposed by @bleyerj in https://gitlab.enpc.fr/navier-fenics/fenics-optim/-/blob/bcd1882f34ffbb9e90fb1df50c184d28daab9162/fenics_optim/custom_quadrature_schemes.py with Quadrature element. Nevertheless, the naive approach does not work, see 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 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
QUAD_DEG = 2
dx_deg = df.dx(metadata={"quadrature_degree":QUAD_DEG})
# Function space with DoFs at quadrature points defined as
# custom monkey patch by @bleyerj
QE_quad = df.FiniteElement("Quadrature", mesh.ufl_cell(),
degree=QUAD_DEG,
quad_scheme="lobatto")
QE = df.FiniteElement(family=fe_space, cell=mesh.ufl_cell(),
degree=fe_degree, quad_scheme="lobatto")
V = df.FunctionSpace(mesh, QE)
# 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.legend()
plt.show()
Dolfin complains:
Exception: Quadrature element must have specified quadrature scheme (lobatto) equal to the integral (default).
What is the proper way to use a custom quadrature scheme?
BTW: Changing to ‘vertex’ quadrature throws the same Exception.
Dolfin 2019.2, Ubuntu LTS 20.04, PPA