TensorElement(...,symmetry=True) is totally broken

[Jan Blechta]

Generated code for element = TensorElement("Discontinuous Lagrange", triangle, 0, symmetry=True) shows that the implementation in FFC is totally wrong.

1. Element has dimension 3 and value size 4. evaluate_basis sets only first three values. Obviously, one of basis functions should contribute to two entries in values.

2. evaluate_dof is wrong as well.

It is a mystery why the element seems to work at leat a bit, see DOLFIN report when both evaluate_basis and evaluate_dof is wrong.

Possibly, this is another candidate for implementing in FIAT.

[Martin Sandve Alnæs]

It's not such a mystery. It's handled correctly in the form compilation but not in the element compilation. Those are separate code paths.

"Possibly, this is another candidate for implementing in FIAT."

I don't see why, the tensor element is basically the same as vector and mixed elements, and the symmetry must be known to the form compiler to detect common subexpressions in forms.

[Jan Blechta]

I don't see why, the tensor element is basically the same as vector and mixed elements, and the symmetry must be known to the form compiler to detect common subexpressions in forms.

The reason is that basis (and dual basis) tabulated by FIAT would already respect symmetry and that would avoid the need for working it around again and again on every place (like now, it's taken into account in form compilation but not in element compilation). But my opinion is not strong here.

[Martin Sandve Alnæs]

For efficient evaluation, ffc should generate code that executes (sum_i u_i phi_i(x)) once and copies the result into output values that by symmetry are equal. That's not fiats job. (Maybe it will be finats job in the future.)

[Jan Blechta]

Sure, but it seems to me that once basis (tabulated by FIAT) is symmetric-tensor-valued (and so dimension is adequately smaller than for non-symmetric tensor element), then there would be no extra evaluation. But maybe I'm wrong here.

[Jan Blechta]

Either this needs a fix or error from evaluate_basis/dof. Otherwise we give to users code which computes rubbish.

[Jan Blechta]

Could be fixable by modifying ffc.representation._generate_foo_offsets to take into account for symmetry.

[Martin Sandve Alnæs]

I tried to fix some symmetry handling in uflacs and hit this bug while testing.

Here's a test snippet that may be useful:

from dolfin import *
from ufl.classes import FloatValue
parameters["form_compiler"]["representation"] = "uflacs"

def test_symmetries():
    mesh = UnitSquareMesh(1, 1)
    element = TensorElement("CG", triangle, 1, symmetry=True)
    V = FunctionSpace(mesh, element)
    f = Function(V)
    c = Constant(((1,2), (2,3)))
    f.interpolate(c)

    fv = f((0.5,0.5))
    value1 = assemble(c**2*dx(mesh))
    value2 = assemble(f**2*dx(mesh))
    print value1  # 18 = 1 + 4 + 4 + 9
    print value2  # 13 = ???
    assert abs(value1 - float(sum(fvc**2 for fvc in [1,2,2,3]))) < 1e-12
    assert abs(value2 - float(sum(fvc**2 for fvc in [1,2,2,3]))) < 1e-12

    # this is wrong, 0 instead of 3 for last component...
    for i in range(len(fv)):
        assert abs(fv[i] - [1.0, 2.0, 2.0, 3.0][i]) < 1e-12

[Martin Sandve Alnæs]

Also look for FIXME comments around physical vs reference value shape/size in representation.py, there are a couple of fixmes pointing to bugs there.

[Martin Sandve Alnæs]

And note that the UFL TensorElement supports generic symmetries, i.e. symmetry[(i,j)] -> (k,l) and not just symmetric-tensor-valued.

[Jan Blechta]

My first attempt was to store ufl.FiniteElement.symmetry() into representation and take this into account in code generation. But that's probably not right. Dealing right with value dimensions and offsets at generation stage seems to need a cleanup, as you say, and that might fix the issue.

UFL symmetries, yes!

Also note how quadrature and tensor deal with symmetry. They use from ufl.permutation import build_component_numbering:

• https://bitbucket.org/fenics-project/ffc/src/32692902865a7623c57c13228caa028f207d8308/ffc/quadrature/quadraturetransformerbase.py?at=master&fileviewer=file-view-default#quadraturetransformerbase.py-880
• https://bitbucket.org/fenics-project/ffc/src/32692902865a7623c57c13228caa028f207d8308/ffc/tensor/monomialtransformation.py?at=master&fileviewer=file-view-default#monomialtransformation.py-469

[Martin Sandve Alnæs]

Good, uflacs also uses build_component_numbering, in uflacs/analysis/modified_terminals.py. This is the code I've tried to clean up a bit: when the transformation to e.g. ReferenceValue(Coefficient(tensor_element)) has been applied, the symmetry dict in the element is no longer relevant, as the element.reference_value_shape() is flattened with the symmetries subtracted from the dimension (see ufl/finiteelement/mixedelement.py).

The element generation should probably also use build_component_numbering so we have one definition of the numbering.

[Jan Blechta]

I believe you're about to fix this, aren't you? Feel free to reassign back to me.

[Jan Blechta]

This is still an issue with the new finite element code.

[Jan Blechta]

It would be good to fix this. It is a bug of type "return junk values without warning".