dolfin team mailing list archive

[Martin Sandve Alnæs <martinal@xxxxxxxxx>]

On 8 June 2011 13:11, Kristian Ølgaard <k.b.oelgaard@xxxxxxxxx> wrote:
> On 8 June 2011 12:11, Martin Sandve Alnæs <martinal@xxxxxxxxx> 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 <martinal@xxxxxxxxx> 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 <martinal@xxxxxxxxx> 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 <745646@xxxxxxxxxxxxxxxxxx>:
>>>>> 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 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], :]
>>>>>
>>>>> _______________________________________________
>>>>> Mailing list: https://launchpad.net/~dolfin
>>>>> Post to     : dolfin@xxxxxxxxxxxxxxxxxxx
>>>>> Unsubscribe : https://launchpad.net/~dolfin
>>>>> More help   : https://help.launchpad.net/ListHelp
>>>>>
>>>>
>>>
>>
>> _______________________________________________
>> Mailing list: https://launchpad.net/~dolfin
>> Post to     : dolfin@xxxxxxxxxxxxxxxxxxx
>> Unsubscribe : https://launchpad.net/~dolfin
>> More help   : https://help.launchpad.net/ListHelp
>>
>