value_shape and quadrature element in subclass expression

Question

Hello,

I'm subclassing a vector-valued Expression, and everything works as expected:

class MyExpr(dolfin.Expression):
    def value_shape(self):
        return (3,)

    def eval(self, value, x):
        value[0] = 0.
        value[1] = 1.
        value[2] = 2.

myexpr = MyExpr()
print myexpr.ufl_shape # -> (3,)

val = [float()]*3    
myexpr.eval(val, [0., 0., 0.])
print val # -> [0.0, 1.0, 2.0]

However, whenever I try to make the expression "pointwise", the evaluation works fine but the shape information seems to be messed up:

fe = dolfin.FiniteElement(
    "Quadrature",
    degree=1)
myexpr2 = MyExpr(element=fe)
print myexpr2.ufl_shape # -> ()

val = [float()]*3
myexpr2.eval(val, [0., 0., 0.])
print val # -> [0.0, 1.0, 2.0]

As a result of this, the first expression can be projected onto a vector function space of dimension 3:

mesh = dolfin.RectangleMesh(
    dolfin.Point(0., 0., 0.),
    dolfin.Point(1., 1., 0.),
    1, 1)
fs = dolfin.VectorFunctionSpace(
    mesh=mesh,
    family="Lagrange",
    degree=1,
    dim=3)
func = dolfin.project(
    myexpr,
    V=fs) # -> works fine

But the second one cannot:

Shape mismatch.
Traceback (most recent call last):
  File "test_shape.py", line 58, in <module>
    func = dolfin.project(
        myexpr2,
        V=fs)
  File "/home/genet/.hashdist/bld/profile/hickh3zjf4k2/lib/python2.7/site-packages/dolfin/fem/projection.py", line 107, in project
    L = ufl.inner(w, v)*dx
  File "/home/genet/.hashdist/bld/profile/hickh3zjf4k2/lib/python2.7/site-packages/ufl/operators.py", line 130, in inner
    return Inner(a, b)
  File "/home/genet/.hashdist/bld/profile/hickh3zjf4k2/lib/python2.7/site-packages/ufl/tensoralgebra.py", line 156, in __new__
    ufl_assert(ash == bsh, "Shape mismatch.")
  File "/home/genet/.hashdist/bld/profile/hickh3zjf4k2/lib/python2.7/site-packages/ufl/assertions.py", line 37, in ufl_assert
    if not condition: error(*message)
  File "/home/genet/.hashdist/bld/profile/hickh3zjf4k2/lib/python2.7/site-packages/ufl/log.py", line 151, in error
    raise self._exception_type(self._format_raw(*message))
ufl.log.UFLException: Shape mismatch.

Is that the intended behavior? If so, why? And, what would be the correct way to use vector-valued pointwise expressions? Thanks!

Answer 1

You need to use FiniteElement with proper shape

fe = VectorElement("Quadrature", cell=mesh.ufl_cell(), degree=1, dim=3)

Comments

It seems that the tensor version of the previous fix does not work:

import dolfin

mesh = dolfin.UnitSquareMesh(1,1)

tfs = dolfin.TensorFunctionSpace(
    mesh=mesh,
    family="Lagrange",
    degree=1)

class MyTensorExpr(dolfin.Expression):
    def value_shape(self):
        return (2,2)

    def eval(self, value, x):
        print "x = " + str(x)
        value[0] = 0.
        value[1] = 1.
        value[2] = 2.
        value[3] = 3.

te_poly = dolfin.TensorElement(
    family="Lagrange",
    cell=mesh.ufl_cell(),
    degree=q_poly,
    quad_scheme="default")

expr_poly = MyTensorExpr(element=te_poly)
print "expr_poly.ufl_shape = " + str(expr_poly.ufl_shape)

func_poly = dolfin.project(
    expr_poly,
    V=tfs,
    form_compiler_parameters={'quadrature_degree':q_poly})
print "func_poly.vector().array() = " + str(func_poly.vector().array())

te_quad = dolfin.TensorElement(
    family="Quadrature",
    cell=mesh.ufl_cell(),
    degree=q_quad,
    quad_scheme="default")

expr_quad = MyTensorExpr(element=te_quad)
print "expr_quad.ufl_shape = " + str(expr_quad.ufl_shape)

func_quad = dolfin.project(
    expr_quad,
    V=tfs,
    form_compiler_parameters={'quadrature_degree':q_quad})
print "func_quad.vector().array() = " + str(func_quad.vector().array())

Here the polynomial version works as expected, but with the quadrature version I'm getting the following Exception:

Quadrature element must have specified quadrature scheme (None) equal to the integral (default).

Am I missing something? Thanks!

---

quad_scheme kwarg does not seem to work as intended for TensorElement. The intefrace is going to change anyway https://bitbucket.org/fenics-project/ufl/issues/78.

---

Thanks! Do you know of any workaround for now?

---

Fiddling with te_quad._quad_scheme (and possibly te_quad.sub_elements()[i]._quad_scheme might help.

---

Awesome! Thanks!