dolfin team mailing list archive

[Kristian Ølgaard <k.b.oelgaard@xxxxxxxxx>]

On 8 June 2011 13:31, Martin Sandve Alnæs <martinal@xxxxxxxxx> wrote:
> 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? :)

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 <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
>>>
>>
>