How to increase the order of the quadrature?

Question

I want to integrate the product of a linear basis function and a characteristic function. For this I want to apply a quadrature rule that considers more points than just the vertices. However, setting parameters["form_compiler_parameters"]["quadrature_degree"] to a higher number, does not change anything. (See the example below)

Since for quadratic basis functions, a higher order scheme is used, my guess is that FEniCS simply sets a quadrature scheme that it considers sufficient.

Is this right? What can I do??

import dolfin


class ContDomain(dolfin.SubDomain):
    """define a subdomain"""
    def __init__(self, ddict):
        dolfin.SubDomain.__init__(self)
        self.minxy = [ddict['xmin'], ddict['ymin']]
        self.maxxy = [ddict['xmax'], ddict['ymax']]

    def inside(self, x, on_boundary):
        return (dolfin.between(x[0], (self.minxy[0], self.maxxy[0]))
                and
                dolfin.between(x[1], (self.minxy[1], self.maxxy[1])))


class CharactFun(dolfin.Expression):
    """ characteristic function of subdomain """
    def __init__(self, subdom):
        self.subdom = subdom
        self.evalcount = 0

    def eval(self, value, x):
        if self.subdom.inside(x, False):
            value[:] = 1
        else:
            value[:] = 0
        self.evalcount += 1

odcoo = dict(xmin=0.25,
             xmax=0.75,
             ymin=0.25,
             ymax=0.75)

mesh = dolfin.UnitSquareMesh(1, 1)
V = dolfin.FunctionSpace(mesh, "CG", 1)
dolfin.plot(mesh)

odom = ContDomain(odcoo)

charfun = CharactFun(odom)
v = dolfin.TestFunction(V)

for qdeg in range(1, 5):
    dolfin.parameters["form_compiler"]["quadrature_degree"] = qdeg

    print 'Quadrature degree: {0}'.\
        format(dolfin.parameters["form_compiler"]["quadrature_degree"])

    vchar = dolfin.assemble(v * charfun * dolfin.dx)

    print 'Values of the integrals: {}'.format(vchar.array())
    print ('Number of evaluations ' +
           'of the charfun: {0}\n').format(charfun.evalcount)
    charfun.evalcount = 0

dolfin.interactive(True)

Gives the output:

~/work/code/optconpy$ python minex_quadraturefenics.py 
Quadrature degree: 1
Calling FFC just-in-time (JIT) compiler, this may take some time.
Values of the integrals: [ 0.  0.  0.  0.]
Number of evaluations of the charfun: 6

Quadrature degree: 2
Calling FFC just-in-time (JIT) compiler, this may take some time.
Values of the integrals: [ 0.  0.  0.  0.]
Number of evaluations of the charfun: 6

Quadrature degree: 3
Values of the integrals: [ 0.  0.  0.  0.]
Number of evaluations of the charfun: 6

Quadrature degree: 4
Values of the integrals: [ 0.  0.  0.  0.]
Number of evaluations of the charfun: 6

Answer 1

I'm sorry for my previous explanation not being fully correct. Expression's element and degree are not used merely for quadrature degree estimation but also for an interpolation of the expression hidden in the DOLFIN-(generated-UFC-code) evaluation chain. Unless set by user, Lagrange of degree 1 is used by default. So your previous interpretation was rather correct.

For your code it means that you call

charfun = CharactFun(odom, degree=expdeg)

with expdeg > 1, possibly with expdeg = qdeg - 1. To get this working, you need to alter your constructor signature to

class CharactFun(dolfin.Expression):
    """ characteristic function of subdomain """
    def __init__(self, subdom, **kwargs):

Note: if you set expression's degree then it is useless to set quadrature degree manually because your integrand has well-defined polynomial degree expdeg + 1 and sufficient quadrature degree should be estimated correctly by UFL algorithms.

Comments

Thank you for this insight. With your modification everything works as I expected.