[firedrake] Index notation with higher-order derivatives

Christian Jacobs c.jacobs10 at imperial.ac.uk
Tue Jun 23 18:57:31 BST 2015


Many thanks Martin - that works for me. I just had to change .dx((j,)*k) to
.dx(*((j,)*k)) in order for the tuple to be used as an arguments list,
otherwise

from firedrake import *
mesh = UnitSquareMesh(10, 10)
tfs = TensorFunctionSpace(mesh, "CG", 1)
vfs = VectorFunctionSpace(mesh, "CG", 1)
w = TestFunction(vfs)
s0 = Function(tfs).interpolate(Expression([[1,2], [3,4]]))
k = 2
dim = s0.ufl_shape[1]
i = Index()
result = sum(w[i] * s0[i,j].dx((j,)*k) for j in range(dim))

fails with

Traceback (most recent call last):
  File "index.py", line 11, in <module>
    result = sum(w[i] * s0[i,j].dx((j,)*k) for j in range(dim))
  File "index.py", line 11, in <genexpr>
    result = sum(w[i] * s0[i,j].dx((j,)*k) for j in range(dim))
  File "/usr/local/lib/python2.7/dist-packages/ufl/exproperators.py", line
418, in _dx
    return d.__getitem__((Ellipsis,) + ii)
  File "/usr/local/lib/python2.7/dist-packages/ufl/exproperators.py", line
367, in _getitem
    all_indices, slice_indices, repeated_indices =
create_slice_indices(component, shape, self.ufl_free_indices)
  File
"/usr/local/lib/python2.7/dist-packages/ufl/index_combination_utils.py",
line 178, in create_slice_indices
    error("Not expecting {0}.".format(ind))
  File "/usr/local/lib/python2.7/dist-packages/ufl/log.py", line 151, in
error
    raise self._exception_type(self._format_raw(*message))
ufl.log.UFLException: Not expecting (0, 0).

Cheers,
Christian

On 23 June 2015 at 10:53, Martin Sandve Alnæs <martinal at simula.no> wrote:

> You're doing something like f.dx(i).dx(i).dx(i), which implies summation
> over the first two i's and leaves the third one free. Try this:
>
> k = some_integer
>
> # First term written explicitly:
> w[0] * sigma[0,0].dx((0,)*k)
>
> # All terms:
> dim = sigma.ufl_shape[1]
> i = Index()
> result = sum(w[i] * sigma[i,j].dx((j,)*k) for j in range(dim))
>
> Martin
>
>
> On 23 June 2015 at 02:05, Christian Jacobs <c.jacobs10 at imperial.ac.uk>
> wrote:
>
>> I'm trying to express this in UFL:
>> http://amcg.ese.ic.ac.uk/~ctj10/images/equations.png
>>
>> For the case of k = 2, I can write
>> F = inner(w[0], grad(grad(s0[0,0])[0])[0] + grad(grad(s0[0,1])[1])[1])*dx
>> + inner(w[1], grad(grad(s0[1,0])[0])[0] + grad(grad(s0[1,1])[1])[1])*dx
>> where w is a TestFunction in a VectorFunctionSpace and s0 (sigma) is a
>> Function in a TensorFunctionSpace.
>>
>> I'd like to generalise this for an arbitrary value of k. I've tried
>>
>> Aij = s0[i,j]
>> k = 2
>> for c in range(0, k):
>>    Aij = Aij.dx(j)
>> A = as_vector(Aij, i)
>> print A
>> F = inner(w, A)*dx
>>
>> For k = 1, { A | A_{i_0} = sum_{i_1} (grad(StressOld[i_0, i_1]))[i_1]  }
>> For k = 2, { A | A_{i_0} = (grad(sum_{i_1} (grad(StressOld[i_0,
>> i_1]))[i_1] ))[i_1] }
>> For k = 3, { A | A_{i_0} = sum_{i_1} (grad((grad(sum_{i_1}
>> (grad(StressOld[i_0, i_1]))[i_1] ))[i_1]))[i_1]
>>
>> The cases k = 1 and k = 3 compile fine, but with k = 2 I'm getting the
>> following error message:
>> ufl.log.UFLException: Can only integrate scalar expressions. The
>> integrand is a tensor expression with value rank 0 and free indices (1,).
>>
>> Using inner(w[i], A) doesn't get rid of the "free indices" issue. Any
>> ideas as to what's going wrong here?
>>
>> The error can be reproduced with:
>>
>> from firedrake import *
>> mesh = UnitSquareMesh(10, 10)
>> tfs = TensorFunctionSpace(mesh, "CG", 1)
>> vfs = VectorFunctionSpace(mesh, "CG", 1)
>> w = TestFunction(vfs)
>> s0 = Function(tfs).interpolate(Expression([[1,2], [3,4]]))
>> Aij = s0[i,j]
>> k = 2
>> for c in range(0, k):
>>    Aij = Aij.dx(j)
>> A = as_vector(Aij, i)
>> print A
>> F = inner(w, A)*dx
>>
>> _______________________________________________
>> firedrake mailing list
>> firedrake at imperial.ac.uk
>> https://mailman.ic.ac.uk/mailman/listinfo/firedrake
>>
>>
>
> _______________________________________________
> firedrake mailing list
> firedrake at imperial.ac.uk
> https://mailman.ic.ac.uk/mailman/listinfo/firedrake
>
>
-------------- next part --------------
HTML attachment scrubbed and removed


More information about the firedrake mailing list