Manually assign tensor field values

Question

I have a Tensor, and I want to replace the values for vertex 2 like this:

vert_2[0,0]=0
vert_2[0,0]=1
vert_2[0,0]=2
vert_2[0,0]=3

is this a correct way to do it?

from dolfin import *

M=UnitSquareMesh(4,4)

TS=TensorFunctionSpace(M, "CG", 1)

v_to_d=TS.dofmap().vertex_to_dof_map(M)
d_to_v=TS.dofmap().dof_to_vertex_map(M)

F=Function(TS)

F=interpolate(Expression((('2*x[0]','x[1]'),('1.0','2.0'))),TS)

ve=F.vector().array()[d_to_v]

n=2

ve[n*4]=0
ve[n*4+1]=1
ve[n*4+2]=2
ve[n*4+3]=3

F.vector().set_local(ve[v_to_d])

With this, I assume that when one does:

ve=F.vector().array()[d_to_v]

he gets back an array that's ordered like this:

vertex 0 component 0,0
vertex 0 component 0,1
vertex 0 component 1,0
vertex 0 component 1,1
...
vertex n component 0,0
vertex n component 0,1
vertex n component 1,0
vertex n component 1,1
...

is it correct?

Thanks.

Answer 1

This seems correct - does it work? Nevertheless it should be much simpler using

F.vector()[d_to_v[4*n:4*(n+1)]] = numpy.array([0.0, 1.0, 2.0, 3.0])

which avoids copying all the DOFs which you don't actually set. But note that there there would be a substantial problems with this in parallel. See Setting Function values in parallel.

Comments

Thanks Jan, that's even better. By the way yes, it works, but I wanted to be sure that the procedure was theoretically correct since the dof / vertex mapping is sometimes obscure to me. In particular I wanted to make sure that the order of the components was always 00 01 10 11 and it din't need any further mapping. Thanks again.

---

make sure that the order of the components was always 00 01 10 11

Well, I'm not sure about this part. I suggest you just test it by something like

u = project(as_matrix(((1.0, 2.0),(3.0, 4.0))), TS)

and then looking on the DOFs by

u.vector()[d_to_v]

---

Ok I did some test on this topic and I'll post here some of them just in case it might be useful for someone else.

1 - Take a tensor field, edit the values for one node, and assemble it again:

from dolfin import *

M=UnitSquareMesh(4,4)

TS=TensorFunctionSpace(M, "CG", 1)

vd=TS.dofmap().vertex_to_dof_map(M)
dv=TS.dofmap().dof_to_vertex_map(M)

F=Function(TS)

F=interpolate(Expression((('2*x[0]','x[1]'),('1.0','2.0'))),TS)

ve=F.vector().array()[dv]

#Change the values for node 2, for instance
n=2

ve[n*4]=0
ve[n*4+1]=1
ve[n*4+2]=2
ve[n*4+3]=3

F.vector().set_local(ve[vd])

#plot one component to see what happened
plot(F[0,0])
interactive()

2 - Manually edit the values for all the vertices of an element:

from dolfin import *
import numpy

M=UnitSquareMesh(4,4)
nv=M.num_vertices()
ce=M.cells()
TS=TensorFunctionSpace(M, "CG", 1)
MET2=Function(TS)

arr=numpy.zeros((nv*4,))
map1=TS.dofmap().vertex_to_dof_map(M)

#set 1234 for all the vertices of element 1
indl=ce[1]
arr[indl*4]=1
arr[indl*4+[1,1,1]]=2
arr[indl*4+[2,2,2]]=3
arr[indl*4+[3,3,3]]=4

#Set the values using the map
MET2.vector().set_local(arr[map1])

a2,b2,c2,d2=MET2.split(deepcopy=True)
a2v=a2.vector().array()

#Print the 1st component ordered as one would expect
print a2v[a2.function_space().dofmap().dof_to_vertex_map(M)]

map2=TS.dofmap().dof_to_vertex_map(M)

#Print the whole tensor as one would expect
D=MET2.vector().array()[map2]

print D

3 - Compute the hessian of a function and correctly plot the values as one would expect to see them:

from dolfin import *
import numpy

M=UnitSquareMesh(4,4)
nv=M.num_vertices()
coo=M.coordinates()
TS=TensorFunctionSpace(M, "CG", 1)

V=FunctionSpace(M,"CG",3)
u=interpolate(Expression('x[0]*x[0]+x[1]*x[1]*x[1]/6.0+x[0]*x[1]'),V)
gradut = project(grad(u), VectorFunctionSpace(M, "CG", 2))
hessut = project(grad(gradut), TS)

arr=numpy.zeros((nv*4,))

mappadv=TS.dofmap().dof_to_vertex_map(M)

arr=hessut.vector().array()[mappadv]

for i in range(0,nv):
    print "vertex %d,coords(%.2f,%.2f)\t"%(i,coo[i,0],coo[i,1]),
    for j in range(0,4):
        print "%.2f"%(arr[i*4+j]),
    print ""