# Custom quadrature rule for Quadrature element

**URL:** <https://fenicsproject.discourse.group/t/custom-quadrature-rule-for-quadrature-element/6597>\
**Category:** Uncategorized\
**Created:** [September 16, 2021, 1:39pm UTC](https://fenicsproject.discourse.group/t/custom-quadrature-rule-for-quadrature-element/6597 "2021-09-16T13:39:47Z")\
**Posts on this page:** 2\
**Page:** 1

<div class="post-metadata">

**Author:** ![petr-henys](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/petr-henys/32/2958_2.png) [@petr-henys](https://fenicsproject.discourse.group/u/petr-henys)\
**Post date:** [September 16, 2021, 1:39pm UTC](https://fenicsproject.discourse.group/t/custom-quadrature-rule-for-quadrature-element/6597/1 "2021-09-16T13:39:47Z")

</div>

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](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:

```auto
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

---

<div class="post-metadata">

**Author:** ![bleyerj](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/bleyerj/32/6020_2.png) [@bleyerj](https://fenicsproject.discourse.group/u/bleyerj)\
**Post date:** [September 17, 2021, 6:58pm UTC](https://fenicsproject.discourse.group/t/custom-quadrature-rule-for-quadrature-element/6597/2 "2021-09-17T18:58:38Z")

</div>

Hi,  
first I think you have to import the custom\_quadrature\_schemes module, then you have to specify the “scheme” to the integral measure:

```auto
dx_deg = df.dx(metadata={"quadrature_degree": QUAD_DEG}, scheme=scheme)

```

I hope this helps.