evaluate_basis error for combination mixed(piola, other) on manifolds

[Martin Sandve Alnæs]

The reason is that reference and physical component offsets are not handled correctly.

This will be fixed in my current work on rewriting the evaluate_basis generation.

[Martin Sandve Alnæs]

Reproduced with this dolfin test

from dolfin import *

def test_evaluate_mixed_piola_mapped_element_on_manifold():
    mesh = UnitCubeMesh(3, 3, 3)
    bmesh = BoundaryMesh(mesh, "exterior")
    bmesh.init_cell_orientations(Expression(("x[0]-0.5", "x[1]-0.5", "x[2]-0.5")))

    U = FunctionSpace(bmesh, "CG", 1)
    V = FunctionSpace(bmesh, "RT", 1)

    # Mix spaces in different order
    UV = U*V
    VU = V*U

    u = Function(U)
    v = Function(V)
    uv = Function(UV)
    vu = Function(VU)

    ue = Expression("x[0]")
    ve = Expression(("x[1]", "x[2]", "x[0]*x[1]"))
    uve = Expression(("x[0]", "x[1]", "x[2]", "x[0]*x[1]"))
    vue = Expression(("x[1]", "x[2]", "x[0]*x[1]", "x[0]"))
    u.interpolate(ue)
    v.interpolate(ve)
    uv.interpolate(uve)
    vu.interpolate(vue)

    u2 = assemble(u**2*dx)
    v2 = assemble(v**2*dx)
    uv2 = assemble(uv**2*dx)
    vu2 = assemble(vu**2*dx)
    assert round(u2 + v2 - uv2, 7) == 0.0
    assert round(u2 + v2 - vu2, 7) == 0.0

    x = (1.0, 1.0, 1.0)
    ux = u(x)
    vx = v(x)
    uvx = uv(x)
    vux = vu(x)
    assert round(ux - uvx[0], 7) == 0.0
    #assert round(ux - vux[3], 7) == 0.0 # This currently fails
    assert all(round(vx[i] - uvx[i+1], 7) == 0.0 for i in range(3))
    #assert all(round(vx[i] - vux[i], 7) == 0.0 for i in range(3)) # This currently fails

test_evaluate_mixed_piola_mapped_element_on_manifold()

[Martin Sandve Alnæs]

Relevant to https://bitbucket.org/fenics-project/ffc/issues/93/replace-code-generation-for-ufc-classes

[Prof Garth Wells]

Removing milestone: 1.7 (automated comment)

[Jan Blechta]

Might be related to #78.

[Martin Sandve Alnæs]

There's a clearly marked FIXME on the physical value shape in representation building. See:

This uses reference value shape from fiat where physical value shape should be used. Needs to get that from ufl subelements: https://bitbucket.org/fenics-project/ffc/src/fd90840b1c078f171e08c066aff0c42aadbf4666/ffc/representation.py?at=master&fileviewer=file-view-default#representation.py-843

Workaround for non-mixed element: https://bitbucket.org/fenics-project/ffc/src/fd90840b1c078f171e08c066aff0c42aadbf4666/ffc/representation.py?at=master&fileviewer=file-view-default#representation.py-852

Fixing this apparently broke tests: https://bitbucket.org/fenics-project/ffc/src/fd90840b1c078f171e08c066aff0c42aadbf4666/ffc/evaluatebasis.py?at=master&fileviewer=file-view-default#evaluatebasis.py-270

[Martin Sandve Alnæs]

This is maybe fixed now?

[Jan Blechta]

Possibly fixed by @garth-wells in #161. Modified Martin's example now works

from dolfin import *

def test_evaluate_mixed_piola_mapped_element_on_manifold():
    mesh = UnitCubeMesh(3, 3, 3)
    bmesh = BoundaryMesh(mesh, "exterior")
    bmesh.init_cell_orientations(Expression(("x[0]-0.5", "x[1]-0.5", "x[2]-0.5"), degree=1))

    U = FunctionSpace(bmesh, "CG", 1)
    V = FunctionSpace(bmesh, "RT", 1)

    # Mix spaces in different order
    UV = FunctionSpace(bmesh, U.ufl_element()*V.ufl_element())
    VU = FunctionSpace(bmesh, V.ufl_element()*U.ufl_element())

    u = Function(U)
    v = Function(V)
    uv = Function(UV)
    vu = Function(VU)

    ue = Expression("x[0]", degree=1)
    ve = Expression(("x[1]", "x[2]", "x[0]*x[1]"), degree=2)
    uve = Expression(("x[0]", "x[1]", "x[2]", "x[0]*x[1]"), degree=2)
    vue = Expression(("x[1]", "x[2]", "x[0]*x[1]", "x[0]"), degree=2)
    u.interpolate(ue)
    v.interpolate(ve)
    uv.interpolate(uve)
    vu.interpolate(vue)

    u2 = assemble(u**2*dx)
    v2 = assemble(v**2*dx)
    uv2 = assemble(uv**2*dx)
    vu2 = assemble(vu**2*dx)
    assert round(u2 + v2 - uv2, 7) == 0.0
    assert round(u2 + v2 - vu2, 7) == 0.0

    x = (1.0, 1.0, 1.0)
    ux = u(x)
    vx = v(x)
    uvx = uv(x)
    vux = vu(x)
    assert round(ux - uvx[0], 7) == 0.0
    assert round(ux - vux[3], 7) == 0.0 # This currently fails
    assert all(round(vx[i] - uvx[i+1], 7) == 0.0 for i in range(3))
    assert all(round(vx[i] - vux[i], 7) == 0.0 for i in range(3)) # This currently fails

test_evaluate_mixed_piola_mapped_element_on_manifold()

[Jan Blechta]

Fixed probably by 080efc82d64b.