On 8 June 2011 13:31, Martin Sandve Alnæs <marti...@simula.no> wrote:
> On 8 June 2011 13:11, Kristian Ølgaard <k.b.oelga...@gmail.com> wrote:
>> On 8 June 2011 12:11, Martin Sandve Alnæs <marti...@simula.no> 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 <marti...@simula.no> 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 <marti...@simula.no> 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 <745...@bugs.launchpad.net>:
>>>>>> 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@lists.launchpad.net
>>>>>> Unsubscribe : https://launchpad.net/~dolfin
>>>>>> More help : https://help.launchpad.net/ListHelp
>>>>>>
>>>>>
>>>>
>>>
>>> _______________________________________________
>>> Mailing list: https://launchpad.net/~dolfin
>>> Post to : dolfin@lists.launchpad.net
>>> Unsubscribe : https://launchpad.net/~dolfin
>>> More help : https://help.launchpad.net/ListHelp
>>>
>>
>
_______________________________________________
Mailing list: https://launchpad.net/~dolfin
Post to : dolfin@lists.launchpad.net
Unsubscribe : https://launchpad.net/~dolfin
More help : https://help.launchpad.net/ListHelp
