[firedrake] Index notation with higher-order derivatives

[Martin Sandve Alnæs]

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
>
>
-------------- next part --------------
HTML attachment scrubbed and removed