Spatial derivative on conditional using vertex quadrature scheme fails

[Matt Covalt]

I posted this over here, but I think this might be a bug instead of user error.

FEniCS fails when evaluating a spatial derivative on a conditional while implementing the vertex quadrature scheme with the following exit message:

ufl.log.UFLException: Illegal component index '(2,)' (value rank 1)for element (value rank 0)​

Here is a MWE:

from fenics import *

# FEniCS primitives
mesh = UnitSquareMesh(10,10)
V = FunctionSpace(mesh, 'CG', 1)
v = TestFunction(V)
u = TrialFunction(V)
sol = Function(V)
sol_prev = interpolate(Constant(0.0), V)  # initial guess

# Nonlinear functions
def q(var):
    return var**2

def q_cond(var):
    return conditional(lt(var, 0), 0, var**2)

quadx = dx(metadata={'quadrature_degree': 1,
                     'quadrature_rule': 'vertex',
                     'representation': 'quadrature'})

#########
# Forms #
#########
# Works: non-conditional w/ default quadrature scheme
a = q(sol_prev)*dot(grad(u), grad(v))*dx
L = q(sol_prev).dx(0)*v*dx
solve(a == L, sol)

# Works: conditional w/ default quadrature scheme
a = q_cond(sol_prev)*dot(grad(u), grad(v))*dx
L = q_cond(sol_prev).dx(0)*v*dx
solve(a == L, sol)

# Works: non-conditional w/ vertex quadrature scheme
a = q(sol_prev)*dot(grad(u), grad(v))*quadx
L = q(sol_prev).dx(0)*v*quadx
solve(a == L, sol)

# Fails: Conditional w/ vertex quadrature scheme
a = q_cond(sol_prev)*dot(grad(u), grad(v))*quadx
L = q_cond(sol_prev).dx(0)*v*quadx
solve(a == L, sol)​

[Jan Blechta]

Yes, that's a bug in quadrature representation which is deprecated and will not be fixed. At the same time there is a bug in uflacs representation: https://bitbucket.org/fenics-project/ffc/issues/145.

It seems it is possible to work around the bug in quadrature representation by propagating the derivative to conditional arguments:

def q_cond(var):
    return conditional(lt(var, 0), 0, (var**2).dx(0))