Piecewise definition of grad (help with syntax)

Question

First, I'm sry for such trivial question (once again) but I cannot google the solution and also didn't find it in the FENICS book (although I didn't read everything). I simplified my original problem to the problem of implementing the grad of vector function.

from dolfin import *

dim = 2

mesh = UnitSquareMesh(5,5)

V = VectorFunctionSpace(mesh, 'DG', 1)
W = TensorFunctionSpace(mesh, 'DG', 1, shape=(dim,dim))

U = Function(V)
M = Function(W)

for i in range(dim-1):
    for j in range(dim-1):
        assign( M[i,j], U[i].dx(j) )

According to FENICS book (pg. 313) the j-th derivative of i-th component of vector function U should be U[i].dx(j). The result I'm getting with this code is that the assign method ends with error 'list of Function expected'.
I also tried the following option

ind = i*dim + j 
assign( M.sub(ind), U[i].dx(j) )

which gave me the same error.

Can you pls help me with the correct syntax? I would also appreciate if you can recomend me some tutorial to syntax in FENICS. How to work with vectors, functions, derivatives - how to acces components of those. Im kind of lost in those and cannot find the answers anywhere.

Answer 1

Hi, you can assign from a Function object or a list of Function objects. The type of U[i].dx(j) in your code is 'ufl.indexed.Indexed - hence the error. Consider instead

from dolfin import *

dim = 2

mesh = UnitSquareMesh(5,5)

S = FunctionSpace(mesh, 'DG', 1)
V = VectorFunctionSpace(mesh, 'DG', 1)
W = TensorFunctionSpace(mesh, 'DG', 1, shape=(dim,dim))

U = interpolate(Expression(('x[0]', 'x[1]')), V)
M = Function(W)

for i in range(dim):
    for j in range(dim):
        print type(U[i].dx(j))
        ind = i*dim + j 
        assign(M.sub(ind), project(U[i].dx(j), S))

for i in range(dim):
    for j in range(dim):
        plot(M[i, j], title='[%d, %d]' % (i, j))
interactive()

Here, a Function (result of projection) is assigned. As for syntax, UFL manual is a good starting point. For components etc, search the forum for component.

Comments

Thank you very much. Nevertheless, is there a faster possibility than using projection? It seeems to me that it is not very fast - in comparison to direct acces to memory. I mean, just in my naive understanding, since U is a vector from space SxS, its components are from space S and more importantly they have to be somehow distinguished in memory where U is saved. There really isnt a way to acces those components directly?