Problem with assemble() with MixedFunctionSpace of symmetric TensorFunctionSpaces

[Jørgen Myre]

from dolfin import *

mesh = Rectangle(0,0,1,1,10,10)
dim = 2 #assume 2D
symm = dict(((i,j), (j,i))
            for i in range(dim) for j in range(dim) if i > j )

T1 = TensorFunctionSpace(mesh, 'CG', 1, symmetry=symm)
T2 = TensorFunctionSpace(mesh, 'CG', 1, symmetry=symm)

TT = T1*T2

R, S = TrialFunctions(TT)
v_R, v_S = TestFunctions(TT)

P = Function(T1)
Q = Function(T2)

Fr = inner(R,v_R)*dx+ inner(P,v_R)*dx

Fs = inner( S, v_S )*dx + inner( Q, v_S )*dx

F = Fr + Fs

A = lhs(F)
B = rhs(F)

a = Matrix()
a = assemble(A, tensor=a)
b = Vector()
b = assemble(B, tensor=b)

I get the following error:

Calling FFC just-in-time (JIT) compiler, this may take some time.
Traceback (most recent call last):
  File "tensortest2.py", line 28, in <module>
    a = assemble(A, tensor=a)
  File "/usr/lib/python2.6/dist-packages/dolfin/fem/assemble.py", line 100, in assemble
    common_cell=common_cell)
  File "/usr/lib/python2.6/dist-packages/dolfin/fem/form.py", line 34, in __init__
    (self._compiled_form, module, self.form_data) = jit(form, form_compiler_parameters, common_cell)
  File "/usr/lib/python2.6/dist-packages/dolfin/compilemodules/jit.py", line 47, in mpi_jit
    return local_jit(*args, **kwargs)
  File "/usr/lib/python2.6/dist-packages/dolfin/compilemodules/jit.py", line 114, in jit
    return jit_compile(form, parameters=p, common_cell=common_cell)
  File "/usr/lib/python2.6/dist-packages/ffc/jitcompiler.py", line 64, in jit
    return jit_form(object, parameters, common_cell)
  File "/usr/lib/python2.6/dist-packages/ffc/jitcompiler.py", line 122, in jit_form
    compile_form(preprocessed_form, prefix=jit_object.signature(), parameters=parameters)
  File "/usr/lib/python2.6/dist-packages/ffc/compiler.py", line 140, in compile_form
    ir = compute_ir(analysis, parameters)
  File "/usr/lib/python2.6/dist-packages/ffc/representation.py", line 66, in compute_ir
    irs = [_compute_integral_ir(f, i, parameters) for (i, f) in enumerate(forms)]
  File "/usr/lib/python2.6/dist-packages/ffc/representation.py", line 186, in _compute_integral_ir
    parameters)
  File "/usr/lib/python2.6/dist-packages/ffc/tensor/tensorrepresentation.py", line 59, in compute_integral_ir
    ir["AK"] = _compute_terms(monomial_form, None, None, domain_type, quadrature_degree)
  File "/usr/lib/python2.6/dist-packages/ffc/tensor/tensorrepresentation.py", line 98, in _compute_terms
    quadrature_degree)
  File "/usr/lib/python2.6/dist-packages/ffc/tensor/referencetensor.py", line 28, in __init__
    self.A0 = integrate(monomial, domain_type, facet0, facet1, quadrature_order)
  File "/usr/lib/python2.6/dist-packages/ffc/tensor/monomialintegration.py", line 50, in integrate
    psis = [_compute_psi(v, table, len(points), domain_type) for v in monomial.arguments]
  File "/usr/lib/python2.6/dist-packages/ffc/tensor/monomialintegration.py", line 169, in _compute_psi
    Psi[component][tuple(dlist)] = etable[dtuple][:, cindex[0].index_range[component], :]

[Jørgen Myre (jorgenmy)]

NOTE: I initially posted this as a Question, if that's wrong of me I can delete that thread.

I am trying to solve for a system of equations with two symmetric tensors as unknowns, and encountering a number of errors with a MixedFunctionSpace consisting of two TensorFunctionSpaces with symmetry when assembling.

I have written a short piece of code which creates the same error:

What's going wrong here?

[Anders Logg (logg)]

On Wed, Mar 30, 2011 at 12:05:27PM -0000, Jørgen Myre wrote:
> NOTE: I initially posted this as a Question, if that's wrong of me I
> can delete that thread.
>
> I am trying to solve for a system of equations with two symmetric
> tensors as unknowns, and encountering a number of errors with a
> MixedFunctionSpace consisting of two TensorFunctionSpaces with symmetry
> when assembling.
>
> I have written a short piece of code which creates the same error:
>
> What's going wrong here?

Thanks for the bug report. Someone will provide an answer. If we're
silent, it just means we're busy with other things. (I'm proofreading
a paper at the moment.)

--
Anders

[Anders Logg (logg)]

On Wed, Mar 30, 2011 at 12:05:27PM -0000, Jørgen Myre wrote:
> NOTE: I initially posted this as a Question, if that's wrong of me I
> can delete that thread.
>
> I am trying to solve for a system of equations with two symmetric
> tensors as unknowns, and encountering a number of errors with a
> MixedFunctionSpace consisting of two TensorFunctionSpaces with symmetry
> when assembling.
>
> I have written a short piece of code which creates the same error:
>
> What's going wrong here?

Looks like a bug in FFC. Try adding this to your code, somewhere at
the top:

parameters["form_compiler"]["representation"] = "quadrature"

--
Anders

[Kristian B. Ølgaard (k.b.oelgaard)]

The following simple form also fails:

T1 = TensorElement('CG', triangle, 1, symmetry=True)
T2 = TensorElement('CG', triangle, 1, symmetry=True)
TT = T1*T2
P, Q = Coefficients(TT)
M = inner(P, Q)*dx
print M

Printing 'M' results in:

{ ([[
  [ (w_0)[0], (w_0)[1] ],
  [ (w_0)[1], (w_0)[3] ]
]]) : ([[
  [ (w_0)[4], (w_0)[5] ],
  [ (w_0)[5], (w_0)[7] ]
]]) } * dx0

which illustrates the problem in FFC for both tensor and quadrature representations.
The component '7' in the ListTensor does not exist in the MixedElement of FFC which, due to symmetry, only contain 6 'unique' subelements (components).
One could argue that UFL should keep track of symmetry when creating the indices of list tensors such that it maps 3->2, 4->3, 5->4 and 7->6.

[Martin Sandve Alnæs (martinal)]

I checked out what UFL does. If you look at expand_indices.py,
you see that this algorithm does handle symmetries. However,
that algorithm is not (and should not be) used by all form compiler
variants. Here's a simplified version of the code:

    def form_argument(self, x):
        if x.shape() == ():
            return x
        else:
            e = x.element()

            # Get component
            c = self.component()

            # Map it through the symmetry mapping
            if isinstance(e, TensorElement):
                s = e.symmetry() or {}
                c = s.get(c, c)

            return x[c]

I want to make two points from this code:

1) The symmetry mapping is currently only available directly from a
TensorElement. Extending symmetries to be a property of all ufl
expressions could be quite a bit of work. If this is a wanted feature,
make a blueprint.

2) When you do have the symmetry mapping at hand, the mapping is
trivial. It should be a quick fix in FFC, just apply the symmetry
mapping like above everywhere you get the component of a form
argument.

Martin

2011/4/5 Kristian B. Ølgaard <email address hidden>:
> The following simple form also fails:
>
> T1 = TensorElement('CG', triangle, 1, symmetry=True)
> T2 = TensorElement('CG', triangle, 1, symmetry=True)
> TT = T1*T2
> P, Q = Coefficients(TT)
> M = inner(P, Q)*dx
> print M
>
> Printing 'M' results in:
>
> { ([[
>  [ (w_0)[0], (w_0)[1] ],
>  [ (w_0)[1], (w_0)[3] ]
> ]]) : ([[
>  [ (w_0)[4], (w_0)[5] ],
>  [ (w_0)[5], (w_0)[7] ]
> ]]) } * dx0
>
> which illustrates the problem in FFC for both tensor and quadrature representations.
> The component '7' in the ListTensor does not exist in the MixedElement of FFC which, due to symmetry, only contain 6 'unique' subelements (components).
> One could argue that UFL should keep track of symmetry when creating the indices of list tensors such that it maps 3->2, 4->3, 5->4 and 7->6.
>
> --
> You received this bug notification because you are a member of DOLFIN
> Team, which is subscribed to DOLFIN.
> https://bugs.launchpad.net/bugs/745646
>
> Title:
>  Problem with assemble() with MixedFunctionSpace of symmetric
>  TensorFunctionSpaces
>
> Status in DOLFIN:
>  New
>
> Bug description:
>  from dolfin import *
>
>  mesh = Rectangle(0,0,1,1,10,10)
>  dim = 2 #assume 2D
>  symm = dict(((i,j), (j,i))
>              for i in range(dim) for j in range(dim) if i > j )
>
>  T1 = TensorFunctionSpace(mesh, 'CG', 1, symmetry=symm)
>  T2 = TensorFunctionSpace(mesh, 'CG', 1, symmetry=symm)
>
>  TT = T1*T2
>
>  R, S = TrialFunctions(TT)
>  v_R, v_S = TestFunctions(TT)
>
>  P = Function(T1)
>  Q = Function(T2)
>
>  Fr = inner(R,v_R)*dx+ inner(P,v_R)*dx
>
>  Fs = inner( S, v_S )*dx + inner( Q, v_S )*dx
>
>  F = Fr + Fs
>
>  A = lhs(F)
>  B = rhs(F)
>
>  a = Matrix()
>  a = assemble(A, tensor=a)
>  b = Vector()
>  b = assemble(B, tensor=b)
>
>  I get the following error:
>
>  Calling FFC just-in-time (JIT) compiler, this may take some time.
>  Traceback (most recent call last):
>    File "tensortest2.py", line 28, in <module>
>      a = assemble(A, tensor=a)
>    File "/usr/lib/python2.6/dist-packages/dolfin/fem/assemble.py", line 100, in...

[Martin Sandve Alnæs (martinal)]

Update: Of course the expand_indices algorithm does not handle your
example where the tensor element is part of a mixed element. I'll
consider adding symmetry as a property of all elements to fix this,
that is much less work than adding it to all expressions.

Martin

On 8 June 2011 10:58, Martin Sandve Alnæs <email address hidden> wrote:
> I checked out what UFL does. If you look at expand_indices.py,
> you see that this algorithm does handle symmetries. However,
> that algorithm is not (and should not be) used by all form compiler
> variants. Here's a simplified version of the code:
>
>    def form_argument(self, x):
>        if x.shape() == ():
>            return x
>        else:
>            e = x.element()
>
>            # Get component
>            c = self.component()
>
>            # Map it through the symmetry mapping
>            if isinstance(e, TensorElement):
>                s = e.symmetry() or {}
>                c = s.get(c, c)
>
>            return x[c]
>
>
> I want to make two points from this code:
>
> 1) The symmetry mapping is currently only available directly from a
> TensorElement. Extending symmetries to be a property of all ufl
> expressions could be quite a bit of work. If this is a wanted feature,
> make a blueprint.
>
> 2) When you do have the symmetry mapping at hand, the mapping is
> trivial. It should be a quick fix in FFC, just apply the symmetry
> mapping like above everywhere you get the component of a form
> argument.
>
> Martin
>
> 2011/4/5 Kristian B. Ølgaard <email address hidden>:
>> The following simple form also fails:
>>
>> T1 = TensorElement('CG', triangle, 1, symmetry=True)
>> T2 = TensorElement('CG', triangle, 1, symmetry=True)
>> TT = T1*T2
>> P, Q = Coefficients(TT)
>> M = inner(P, Q)*dx
>> print M
>>
>> Printing 'M' results in:
>>
>> { ([[
>>  [ (w_0)[0], (w_0)[1] ],
>>  [ (w_0)[1], (w_0)[3] ]
>> ]]) : ([[
>>  [ (w_0)[4], (w_0)[5] ],
>>  [ (w_0)[5], (w_0)[7] ]
>> ]]) } * dx0
>>
>> which illustrates the problem in FFC for both tensor and quadrature representations.
>> The component '7' in the ListTensor does not exist in the MixedElement of FFC which, due to symmetry, only contain 6 'unique' subelements (components).
>> One could argue that UFL should keep track of symmetry when creating the indices of list tensors such that it maps 3->2, 4->3, 5->4 and 7->6.
>>
>> --
>> You received this bug notification because you are a member of DOLFIN
>> Team, which is subscribed to DOLFIN.
>> https://bugs.launchpad.net/bugs/745646
>>
>> Title:
>>  Problem with assemble() with MixedFunctionSpace of symmetric
>>  TensorFunctionSpaces
>>
>> Status in DOLFIN:
>>  New
>>
>> Bug description:
>>  from dolfin import *
>>
>>  mesh = Rectangle(0,0,1,1,10,10)
>>  dim = 2 #assume 2D
>>  symm = dict(((i,j), (j,i))
>>              for i in range(dim) for j in range(dim) if i > j )
>>
>>  T1 = TensorFunctionSpace(mesh, 'CG', 1, symmetry=symm)
>>  T2 = TensorFunctionSpace(mesh, 'CG', 1, symmetry=symm)
>>
>>  TT = T1*T2
>>
>>  R, S = TrialFunctions(TT)
>>  v_R, v_S = TestFunctions(TT)
>>
>>  P = Function(T1)
>>  Q = Function(T2)
>>
>>  Fr = inner(R,v_R)*dx+ inner(P,v_R)*dx
>>
>>  Fs = inner( S, v_S )*dx + inner( Q, v_S ...

[Martin Sandve Alnæs (martinal)]

Done and checked in. If someone updates FFC to support this, we can
hopefully close this bug.

Martin

On 8 June 2011 11:17, Martin Sandve Alnæs <email address hidden> wrote:
> Update: Of course the expand_indices algorithm does not handle your
> example where the tensor element is part of a mixed element. I'll
> consider adding symmetry as a property of all elements to fix this,
> that is much less work than adding it to all expressions.
>
> Martin
>
> On 8 June 2011 10:58, Martin Sandve Alnæs <email address hidden> wrote:
>> I checked out what UFL does. If you look at expand_indices.py,
>> you see that this algorithm does handle symmetries. However,
>> that algorithm is not (and should not be) used by all form compiler
>> variants. Here's a simplified version of the code:
>>
>>    def form_argument(self, x):
>>        if x.shape() == ():
>>            return x
>>        else:
>>            e = x.element()
>>
>>            # Get component
>>            c = self.component()
>>
>>            # Map it through the symmetry mapping
>>            if isinstance(e, TensorElement):
>>                s = e.symmetry() or {}
>>                c = s.get(c, c)
>>
>>            return x[c]
>>
>>
>> I want to make two points from this code:
>>
>> 1) The symmetry mapping is currently only available directly from a
>> TensorElement. Extending symmetries to be a property of all ufl
>> expressions could be quite a bit of work. If this is a wanted feature,
>> make a blueprint.
>>
>> 2) When you do have the symmetry mapping at hand, the mapping is
>> trivial. It should be a quick fix in FFC, just apply the symmetry
>> mapping like above everywhere you get the component of a form
>> argument.
>>
>> Martin
>>
>> 2011/4/5 Kristian B. Ølgaard <email address hidden>:
>>> The following simple form also fails:
>>>
>>> T1 = TensorElement('CG', triangle, 1, symmetry=True)
>>> T2 = TensorElement('CG', triangle, 1, symmetry=True)
>>> TT = T1*T2
>>> P, Q = Coefficients(TT)
>>> M = inner(P, Q)*dx
>>> print M
>>>
>>> Printing 'M' results in:
>>>
>>> { ([[
>>>  [ (w_0)[0], (w_0)[1] ],
>>>  [ (w_0)[1], (w_0)[3] ]
>>> ]]) : ([[
>>>  [ (w_0)[4], (w_0)[5] ],
>>>  [ (w_0)[5], (w_0)[7] ]
>>> ]]) } * dx0
>>>
>>> which illustrates the problem in FFC for both tensor and quadrature representations.
>>> The component '7' in the ListTensor does not exist in the MixedElement of FFC which, due to symmetry, only contain 6 'unique' subelements (components).
>>> One could argue that UFL should keep track of symmetry when creating the indices of list tensors such that it maps 3->2, 4->3, 5->4 and 7->6.
>>>
>>> --
>>> You received this bug notification because you are a member of DOLFIN
>>> Team, which is subscribed to DOLFIN.
>>> https://bugs.launchpad.net/bugs/745646
>>>
>>> Title:
>>>  Problem with assemble() with MixedFunctionSpace of symmetric
>>>  TensorFunctionSpaces
>>>
>>> Status in DOLFIN:
>>>  New
>>>
>>> Bug description:
>>>  from dolfin import *
>>>
>>>  mesh = Rectangle(0,0,1,1,10,10)
>>>  dim = 2 #assume 2D
>>>  symm = dict(((i,j), (j,i))
>>>              for i in range(dim) for j in range(dim) if i > j )
>>>
>>>  T1 = TensorFunctionSpace(mesh, 'CG', 1, symmetry=symm)
>>>  T2 = T...

[Kristian B. Ølgaard (k.b.oelgaard)]

On 8 June 2011 12:11, Martin Sandve Alnæs <email address hidden> wrote:
> Done and checked in. If someone updates FFC to support this, we can
> hopefully close this bug.

I'm not sure this is enough to handle the bug. If you look at the
output of printing M in the example code I posted you'll see that the
list tensor contains component '7'. This is what you'll get from
calling self.component(), but the
TT.symmetry() only contains:
{(2,): (1,), (6,): (5,)}

Is there some function we need to call first to map the component '7'
--> '6', before looking at symmetries to map '6' --> '5'?
Doing so will get us into trouble with mapping '3' --> '2' since
symmetry will map that to '1'.
The TT element has 2 x 3 sub elements due to symmetry.

Kristian

> Martin
>
>
> On 8 June 2011 11:17, Martin Sandve Alnæs <email address hidden> wrote:
>> Update: Of course the expand_indices algorithm does not handle your
>> example where the tensor element is part of a mixed element. I'll
>> consider adding symmetry as a property of all elements to fix this,
>> that is much less work than adding it to all expressions.
>>
>> Martin
>>
>> On 8 June 2011 10:58, Martin Sandve Alnæs <email address hidden> wrote:
>>> I checked out what UFL does. If you look at expand_indices.py,
>>> you see that this algorithm does handle symmetries. However,
>>> that algorithm is not (and should not be) used by all form compiler
>>> variants. Here's a simplified version of the code:
>>>
>>>    def form_argument(self, x):
>>>        if x.shape() == ():
>>>            return x
>>>        else:
>>>            e = x.element()
>>>
>>>            # Get component
>>>            c = self.component()
>>>
>>>            # Map it through the symmetry mapping
>>>            if isinstance(e, TensorElement):
>>>                s = e.symmetry() or {}
>>>                c = s.get(c, c)
>>>
>>>            return x[c]
>>>
>>>
>>> I want to make two points from this code:
>>>
>>> 1) The symmetry mapping is currently only available directly from a
>>> TensorElement. Extending symmetries to be a property of all ufl
>>> expressions could be quite a bit of work. If this is a wanted feature,
>>> make a blueprint.
>>>
>>> 2) When you do have the symmetry mapping at hand, the mapping is
>>> trivial. It should be a quick fix in FFC, just apply the symmetry
>>> mapping like above everywhere you get the component of a form
>>> argument.
>>>
>>> Martin
>>>
>>> 2011/4/5 Kristian B. Ølgaard <email address hidden>:
>>>> The following simple form also fails:
>>>>
>>>> T1 = TensorElement('CG', triangle, 1, symmetry=True)
>>>> T2 = TensorElement('CG', triangle, 1, symmetry=True)
>>>> TT = T1*T2
>>>> P, Q = Coefficients(TT)
>>>> M = inner(P, Q)*dx
>>>> print M
>>>>
>>>> Printing 'M' results in:
>>>>
>>>> { ([[
>>>>  [ (w_0)[0], (w_0)[1] ],
>>>>  [ (w_0)[1], (w_0)[3] ]
>>>> ]]) : ([[
>>>>  [ (w_0)[4], (w_0)[5] ],
>>>>  [ (w_0)[5], (w_0)[7] ]
>>>> ]]) } * dx0
>>>>
>>>> which illustrates the problem in FFC for both tensor and quadrature representations.
>>>> The component '7' in the ListTensor does not exist in the MixedElement of FFC which, due to symmetry, only contain 6 'unique' subelements (components).
>>>> One could argue t...

[Martin Sandve Alnæs (martinal)]

On 8 June 2011 13:11, Kristian Ølgaard <email address hidden> wrote:
> On 8 June 2011 12:11, Martin Sandve Alnæs <email address hidden> wrote:
>> Done and checked in. If someone updates FFC to support this, we can
>> hopefully close this bug.
>
> I'm not sure this is enough to handle the bug. If you look at the
> output of printing M in the example code I posted you'll see that the
> list tensor contains component '7'. This is what you'll get from
> calling self.component(), but the
> TT.symmetry() only contains:
> {(2,): (1,), (6,): (5,)}
>
> Is there some function we need to call first to map the component '7'
> --> '6', before looking at symmetries to map '6' --> '5'?
> Doing so will get us into trouble with mapping '3' --> '2' since
> symmetry will map that to '1'.
> The TT element has 2 x 3 sub elements due to symmetry.

The 7 is an index into the value index space of the coefficient and is
correct. It has nothing (directly) to do with subelement indexing. I
think you're assuming a closer relation between them than there is?
Let me try to clear it up...

The value index space is contiguous from the point of view of UFL
expressions, but has holes when symmetries are considered. The
noncontiguous value index space will then need to be mapped to a
contiguous subelement space by associating each value index that is
not in the symmetry mapping with a subelement index.

1) We have a component/value index
2) We map that value index to another value index using a symmetry
mapping (e.g. 6->5 and 7->7 in your example)
3) We map from the noncontiguous value index space to the contiguous
subelement index space

Clear as mud? :)

UFL handles (2) for you only when you apply expand_indices.

FFC will have to handle (3) when generating code, it doesn't touch
anything UFL needs to know about. I'll see if I can whip up a quick
utility function for it though.

Martin

>
> Kristian
>
>> Martin
>>
>>
>> On 8 June 2011 11:17, Martin Sandve Alnæs <email address hidden> wrote:
>>> Update: Of course the expand_indices algorithm does not handle your
>>> example where the tensor element is part of a mixed element. I'll
>>> consider adding symmetry as a property of all elements to fix this,
>>> that is much less work than adding it to all expressions.
>>>
>>> Martin
>>>
>>> On 8 June 2011 10:58, Martin Sandve Alnæs <email address hidden> wrote:
>>>> I checked out what UFL does. If you look at expand_indices.py,
>>>> you see that this algorithm does handle symmetries. However,
>>>> that algorithm is not (and should not be) used by all form compiler
>>>> variants. Here's a simplified version of the code:
>>>>
>>>>    def form_argument(self, x):
>>>>        if x.shape() == ():
>>>>            return x
>>>>        else:
>>>>            e = x.element()
>>>>
>>>>            # Get component
>>>>            c = self.component()
>>>>
>>>>            # Map it through the symmetry mapping
>>>>            if isinstance(e, TensorElement):
>>>>                s = e.symmetry() or {}
>>>>                c = s.get(c, c)
>>>>
>>>>            return x[c]
>>>>
>>>>
>>>> I want to make two points from this code:
>>>>
>>>> 1) The symmetry mapping is currently only available directly from a
...

[Kristian B. Ølgaard (k.b.oelgaard)]

On 8 June 2011 13:31, Martin Sandve Alnæs <email address hidden> wrote:
> On 8 June 2011 13:11, Kristian Ølgaard <email address hidden> wrote:
>> On 8 June 2011 12:11, Martin Sandve Alnæs <email address hidden> wrote:
>>> Done and checked in. If someone updates FFC to support this, we can
>>> hopefully close this bug.
>>
>> I'm not sure this is enough to handle the bug. If you look at the
>> output of printing M in the example code I posted you'll see that the
>> list tensor contains component '7'. This is what you'll get from
>> calling self.component(), but the
>> TT.symmetry() only contains:
>> {(2,): (1,), (6,): (5,)}
>>
>> Is there some function we need to call first to map the component '7'
>> --> '6', before looking at symmetries to map '6' --> '5'?
>> Doing so will get us into trouble with mapping '3' --> '2' since
>> symmetry will map that to '1'.
>> The TT element has 2 x 3 sub elements due to symmetry.
>
> The 7 is an index into the value index space of the coefficient and is
> correct. It has nothing (directly) to do with subelement indexing. I
> think you're assuming a closer relation between them than there is?
> Let me try to clear it up...
>
> The value index space is contiguous from the point of view of UFL
> expressions, but has holes when symmetries are considered. The
> noncontiguous value index space will then need to be mapped to a
> contiguous subelement space by associating each value index that is
> not in the symmetry mapping with a subelement index.
>
> 1) We have a component/value index
> 2) We map that value index to another value index using a symmetry
> mapping (e.g. 6->5 and 7->7 in your example)
> 3) We map from the noncontiguous value index space to the contiguous
> subelement index space
>
> Clear as mud? :)

Yes, but since we only deal with the (sub)elements that are actually
present in FFC, it's really inconvenient that we can't get a mapping
from the component to the subelement.
I somehow suspected the FiniteElement.extract_component() to do this,
but it turns out not to be the case.

> UFL handles (2) for you only when you apply expand_indices.
>
> FFC will have to handle (3) when generating code, it doesn't touch
> anything UFL needs to know about. I'll see if I can whip up a quick
> utility function for it though.

That would really be nice.

Kristian

> Martin
>
>>
>> Kristian
>>
>>> Martin
>>>
>>>
>>> On 8 June 2011 11:17, Martin Sandve Alnæs <email address hidden> wrote:
>>>> Update: Of course the expand_indices algorithm does not handle your
>>>> example where the tensor element is part of a mixed element. I'll
>>>> consider adding symmetry as a property of all elements to fix this,
>>>> that is much less work than adding it to all expressions.
>>>>
>>>> Martin
>>>>
>>>> On 8 June 2011 10:58, Martin Sandve Alnæs <email address hidden> wrote:
>>>>> I checked out what UFL does. If you look at expand_indices.py,
>>>>> you see that this algorithm does handle symmetries. However,
>>>>> that algorithm is not (and should not be) used by all form compiler
>>>>> variants. Here's a simplified version of the code:
>>>>>
>>>>>    def form_argument(self, x):
>>>>>        if x.shape() == ():
>>>>>            return x
...

[Martin Sandve Alnæs (martinal)]

On 8 June 2011 13:46, Kristian Ølgaard <email address hidden> wrote:
> On 8 June 2011 13:31, Martin Sandve Alnæs <email address hidden> wrote:
>> On 8 June 2011 13:11, Kristian Ølgaard <email address hidden> wrote:
>>> On 8 June 2011 12:11, Martin Sandve Alnæs <email address hidden> wrote:
>>>> Done and checked in. If someone updates FFC to support this, we can
>>>> hopefully close this bug.
>>>
>>> I'm not sure this is enough to handle the bug. If you look at the
>>> output of printing M in the example code I posted you'll see that the
>>> list tensor contains component '7'. This is what you'll get from
>>> calling self.component(), but the
>>> TT.symmetry() only contains:
>>> {(2,): (1,), (6,): (5,)}
>>>
>>> Is there some function we need to call first to map the component '7'
>>> --> '6', before looking at symmetries to map '6' --> '5'?
>>> Doing so will get us into trouble with mapping '3' --> '2' since
>>> symmetry will map that to '1'.
>>> The TT element has 2 x 3 sub elements due to symmetry.
>>
>> The 7 is an index into the value index space of the coefficient and is
>> correct. It has nothing (directly) to do with subelement indexing. I
>> think you're assuming a closer relation between them than there is?
>> Let me try to clear it up...
>>
>> The value index space is contiguous from the point of view of UFL
>> expressions, but has holes when symmetries are considered. The
>> noncontiguous value index space will then need to be mapped to a
>> contiguous subelement space by associating each value index that is
>> not in the symmetry mapping with a subelement index.
>>
>> 1) We have a component/value index
>> 2) We map that value index to another value index using a symmetry
>> mapping (e.g. 6->5 and 7->7 in your example)
>> 3) We map from the noncontiguous value index space to the contiguous
>> subelement index space
>>
>> Clear as mud? :)
>
> Yes, but since we only deal with the (sub)elements that are actually
> present in FFC, it's really inconvenient that we can't get a mapping
> from the component to the subelement.
> I somehow suspected the FiniteElement.extract_component() to do this,
> but it turns out not to be the case.
>
>> UFL handles (2) for you only when you apply expand_indices.
>>
>> FFC will have to handle (3) when generating code, it doesn't touch
>> anything UFL needs to know about. I'll see if I can whip up a quick
>> utility function for it though.
>
> That would really be nice.

Done :) The latest patch contains code and test showing usage.

But maybe it should be integrated into extract_component, I'll take a
look at that now that I'm into this.

Martin

[Martin Sandve Alnæs (martinal)]

On 8 June 2011 13:55, Martin Sandve Alnæs <email address hidden> wrote:
> On 8 June 2011 13:46, Kristian Ølgaard <email address hidden> wrote:
>> On 8 June 2011 13:31, Martin Sandve Alnæs <email address hidden> wrote:
>>> On 8 June 2011 13:11, Kristian Ølgaard <email address hidden> wrote:
>>>> On 8 June 2011 12:11, Martin Sandve Alnæs <email address hidden> wrote:
>>>>> Done and checked in. If someone updates FFC to support this, we can
>>>>> hopefully close this bug.
>>>>
>>>> I'm not sure this is enough to handle the bug. If you look at the
>>>> output of printing M in the example code I posted you'll see that the
>>>> list tensor contains component '7'. This is what you'll get from
>>>> calling self.component(), but the
>>>> TT.symmetry() only contains:
>>>> {(2,): (1,), (6,): (5,)}
>>>>
>>>> Is there some function we need to call first to map the component '7'
>>>> --> '6', before looking at symmetries to map '6' --> '5'?
>>>> Doing so will get us into trouble with mapping '3' --> '2' since
>>>> symmetry will map that to '1'.
>>>> The TT element has 2 x 3 sub elements due to symmetry.
>>>
>>> The 7 is an index into the value index space of the coefficient and is
>>> correct. It has nothing (directly) to do with subelement indexing. I
>>> think you're assuming a closer relation between them than there is?
>>> Let me try to clear it up...
>>>
>>> The value index space is contiguous from the point of view of UFL
>>> expressions, but has holes when symmetries are considered. The
>>> noncontiguous value index space will then need to be mapped to a
>>> contiguous subelement space by associating each value index that is
>>> not in the symmetry mapping with a subelement index.
>>>
>>> 1) We have a component/value index
>>> 2) We map that value index to another value index using a symmetry
>>> mapping (e.g. 6->5 and 7->7 in your example)
>>> 3) We map from the noncontiguous value index space to the contiguous
>>> subelement index space
>>>
>>> Clear as mud? :)
>>
>> Yes, but since we only deal with the (sub)elements that are actually
>> present in FFC, it's really inconvenient that we can't get a mapping
>> from the component to the subelement.
>> I somehow suspected the FiniteElement.extract_component() to do this,
>> but it turns out not to be the case.
>>
>>> UFL handles (2) for you only when you apply expand_indices.
>>>
>>> FFC will have to handle (3) when generating code, it doesn't touch
>>> anything UFL needs to know about. I'll see if I can whip up a quick
>>> utility function for it though.
>>
>> That would really be nice.
>
> Done :) The latest patch contains code and test showing usage.
>
> But maybe it should be integrated into extract_component, I'll take a
> look at that now that I'm into this.

I applied the symmetry mapping inside extract_component for tensor elements,
that way you shouldn't have to do the symmetry mapping in addition to
calling extract_component. The numbering utility I checked in could still
come in handy though, so I'll let it stay.

Note that when I pushed now I pushed the float formatting fix as well,
so some FFC references will probably need regeneration again :)

Martin

[Kristian B. Ølgaard (k.b.oelgaard)]

On 8 June 2011 14:08, Martin Sandve Alnæs <email address hidden> wrote:
> On 8 June 2011 13:55, Martin Sandve Alnæs <email address hidden> wrote:
>> On 8 June 2011 13:46, Kristian Ølgaard <email address hidden> wrote:
>>> On 8 June 2011 13:31, Martin Sandve Alnæs <email address hidden> wrote:
>>>> On 8 June 2011 13:11, Kristian Ølgaard <email address hidden> wrote:
>>>>> On 8 June 2011 12:11, Martin Sandve Alnæs <email address hidden> wrote:
>>>>>> Done and checked in. If someone updates FFC to support this, we can
>>>>>> hopefully close this bug.
>>>>>
>>>>> I'm not sure this is enough to handle the bug. If you look at the
>>>>> output of printing M in the example code I posted you'll see that the
>>>>> list tensor contains component '7'. This is what you'll get from
>>>>> calling self.component(), but the
>>>>> TT.symmetry() only contains:
>>>>> {(2,): (1,), (6,): (5,)}
>>>>>
>>>>> Is there some function we need to call first to map the component '7'
>>>>> --> '6', before looking at symmetries to map '6' --> '5'?
>>>>> Doing so will get us into trouble with mapping '3' --> '2' since
>>>>> symmetry will map that to '1'.
>>>>> The TT element has 2 x 3 sub elements due to symmetry.
>>>>
>>>> The 7 is an index into the value index space of the coefficient and is
>>>> correct. It has nothing (directly) to do with subelement indexing. I
>>>> think you're assuming a closer relation between them than there is?
>>>> Let me try to clear it up...
>>>>
>>>> The value index space is contiguous from the point of view of UFL
>>>> expressions, but has holes when symmetries are considered. The
>>>> noncontiguous value index space will then need to be mapped to a
>>>> contiguous subelement space by associating each value index that is
>>>> not in the symmetry mapping with a subelement index.
>>>>
>>>> 1) We have a component/value index
>>>> 2) We map that value index to another value index using a symmetry
>>>> mapping (e.g. 6->5 and 7->7 in your example)
>>>> 3) We map from the noncontiguous value index space to the contiguous
>>>> subelement index space
>>>>
>>>> Clear as mud? :)
>>>
>>> Yes, but since we only deal with the (sub)elements that are actually
>>> present in FFC, it's really inconvenient that we can't get a mapping
>>> from the component to the subelement.
>>> I somehow suspected the FiniteElement.extract_component() to do this,
>>> but it turns out not to be the case.
>>>
>>>> UFL handles (2) for you only when you apply expand_indices.
>>>>
>>>> FFC will have to handle (3) when generating code, it doesn't touch
>>>> anything UFL needs to know about. I'll see if I can whip up a quick
>>>> utility function for it though.
>>>
>>> That would really be nice.
>>
>> Done :) The latest patch contains code and test showing usage.
>>
>> But maybe it should be integrated into extract_component, I'll take a
>> look at that now that I'm into this.
>
> I applied the symmetry mapping inside extract_component for tensor elements,
> that way you shouldn't have to do the symmetry mapping in addition to
> calling extract_component. The numbering utility I checked in could still
> come in handy though, so I'll let it stay.

OK, will you add the numbering utility func...

[Martin Sandve Alnæs (martinal)]

On 8 June 2011 14:18, Kristian Ølgaard <email address hidden> wrote:
> On 8 June 2011 14:08, Martin Sandve Alnæs <email address hidden> wrote:
>> On 8 June 2011 13:55, Martin Sandve Alnæs <email address hidden> wrote:
>>> On 8 June 2011 13:46, Kristian Ølgaard <email address hidden> wrote:
>>>> On 8 June 2011 13:31, Martin Sandve Alnæs <email address hidden> wrote:
>>>>> On 8 June 2011 13:11, Kristian Ølgaard <email address hidden> wrote:
>>>>>> On 8 June 2011 12:11, Martin Sandve Alnæs <email address hidden> wrote:
>>>>>>> Done and checked in. If someone updates FFC to support this, we can
>>>>>>> hopefully close this bug.
>>>>>>
>>>>>> I'm not sure this is enough to handle the bug. If you look at the
>>>>>> output of printing M in the example code I posted you'll see that the
>>>>>> list tensor contains component '7'. This is what you'll get from
>>>>>> calling self.component(), but the
>>>>>> TT.symmetry() only contains:
>>>>>> {(2,): (1,), (6,): (5,)}
>>>>>>
>>>>>> Is there some function we need to call first to map the component '7'
>>>>>> --> '6', before looking at symmetries to map '6' --> '5'?
>>>>>> Doing so will get us into trouble with mapping '3' --> '2' since
>>>>>> symmetry will map that to '1'.
>>>>>> The TT element has 2 x 3 sub elements due to symmetry.
>>>>>
>>>>> The 7 is an index into the value index space of the coefficient and is
>>>>> correct. It has nothing (directly) to do with subelement indexing. I
>>>>> think you're assuming a closer relation between them than there is?
>>>>> Let me try to clear it up...
>>>>>
>>>>> The value index space is contiguous from the point of view of UFL
>>>>> expressions, but has holes when symmetries are considered. The
>>>>> noncontiguous value index space will then need to be mapped to a
>>>>> contiguous subelement space by associating each value index that is
>>>>> not in the symmetry mapping with a subelement index.
>>>>>
>>>>> 1) We have a component/value index
>>>>> 2) We map that value index to another value index using a symmetry
>>>>> mapping (e.g. 6->5 and 7->7 in your example)
>>>>> 3) We map from the noncontiguous value index space to the contiguous
>>>>> subelement index space
>>>>>
>>>>> Clear as mud? :)
>>>>
>>>> Yes, but since we only deal with the (sub)elements that are actually
>>>> present in FFC, it's really inconvenient that we can't get a mapping
>>>> from the component to the subelement.
>>>> I somehow suspected the FiniteElement.extract_component() to do this,
>>>> but it turns out not to be the case.
>>>>
>>>>> UFL handles (2) for you only when you apply expand_indices.
>>>>>
>>>>> FFC will have to handle (3) when generating code, it doesn't touch
>>>>> anything UFL needs to know about. I'll see if I can whip up a quick
>>>>> utility function for it though.
>>>>
>>>> That would really be nice.
>>>
>>> Done :) The latest patch contains code and test showing usage.
>>>
>>> But maybe it should be integrated into extract_component, I'll take a
>>> look at that now that I'm into this.
>>
>> I applied the symmetry mapping inside extract_component for tensor elements,
>> that way you shouldn't have to do the symmetry mapping in addition to
>> calling extract_component. T...

[Kristian B. Ølgaard (k.b.oelgaard)]

On 8 June 2011 14:29, Martin Sandve Alnæs <email address hidden> wrote:
> On 8 June 2011 14:18, Kristian Ølgaard <email address hidden> wrote:
>> On 8 June 2011 14:08, Martin Sandve Alnæs <email address hidden> wrote:
>>> On 8 June 2011 13:55, Martin Sandve Alnæs <email address hidden> wrote:
>>>> On 8 June 2011 13:46, Kristian Ølgaard <email address hidden> wrote:
>>>>> On 8 June 2011 13:31, Martin Sandve Alnæs <email address hidden> wrote:
>>>>>> On 8 June 2011 13:11, Kristian Ølgaard <email address hidden> wrote:
>>>>>>> On 8 June 2011 12:11, Martin Sandve Alnæs <email address hidden> wrote:
>>>>>>>> Done and checked in. If someone updates FFC to support this, we can
>>>>>>>> hopefully close this bug.
>>>>>>>
>>>>>>> I'm not sure this is enough to handle the bug. If you look at the
>>>>>>> output of printing M in the example code I posted you'll see that the
>>>>>>> list tensor contains component '7'. This is what you'll get from
>>>>>>> calling self.component(), but the
>>>>>>> TT.symmetry() only contains:
>>>>>>> {(2,): (1,), (6,): (5,)}
>>>>>>>
>>>>>>> Is there some function we need to call first to map the component '7'
>>>>>>> --> '6', before looking at symmetries to map '6' --> '5'?
>>>>>>> Doing so will get us into trouble with mapping '3' --> '2' since
>>>>>>> symmetry will map that to '1'.
>>>>>>> The TT element has 2 x 3 sub elements due to symmetry.
>>>>>>
>>>>>> The 7 is an index into the value index space of the coefficient and is
>>>>>> correct. It has nothing (directly) to do with subelement indexing. I
>>>>>> think you're assuming a closer relation between them than there is?
>>>>>> Let me try to clear it up...
>>>>>>
>>>>>> The value index space is contiguous from the point of view of UFL
>>>>>> expressions, but has holes when symmetries are considered. The
>>>>>> noncontiguous value index space will then need to be mapped to a
>>>>>> contiguous subelement space by associating each value index that is
>>>>>> not in the symmetry mapping with a subelement index.
>>>>>>
>>>>>> 1) We have a component/value index
>>>>>> 2) We map that value index to another value index using a symmetry
>>>>>> mapping (e.g. 6->5 and 7->7 in your example)
>>>>>> 3) We map from the noncontiguous value index space to the contiguous
>>>>>> subelement index space
>>>>>>
>>>>>> Clear as mud? :)
>>>>>
>>>>> Yes, but since we only deal with the (sub)elements that are actually
>>>>> present in FFC, it's really inconvenient that we can't get a mapping
>>>>> from the component to the subelement.
>>>>> I somehow suspected the FiniteElement.extract_component() to do this,
>>>>> but it turns out not to be the case.
>>>>>
>>>>>> UFL handles (2) for you only when you apply expand_indices.
>>>>>>
>>>>>> FFC will have to handle (3) when generating code, it doesn't touch
>>>>>> anything UFL needs to know about. I'll see if I can whip up a quick
>>>>>> utility function for it though.
>>>>>
>>>>> That would really be nice.
>>>>
>>>> Done :) The latest patch contains code and test showing usage.
>>>>
>>>> But maybe it should be integrated into extract_component, I'll take a
>>>> look at that now that I'm into this.
>>>
>>> I applied the symmetry mapping inside extract_co...

[Kristian B. Ølgaard (k.b.oelgaard)]

A fix has been committed such that the following forms:

T1 = TensorElement("CG", triangle, 1, symmetry=True)
T2 = TensorElement("CG", triangle, 1, symmetry=True)
TT = T1*T2
P, Q = Coefficients(TT)
M = inner(P, Q)*dx

and

T1 = VectorElement("CG", triangle, 1, 3)
T2 = VectorElement("CG", triangle, 1, 3)
TT = T1*T2
P, Q = Coefficients(TT)
def s(v):
    return as_matrix([[v[0], v[1]], [v[1], v[2]]])
M = inner(s(P), s(Q))*dx

result in the exact same output for tabulate_tensor() for both quadrature and tensor representations