How to interpolate data at vertices of (3D) cells?

Question

I am trying to get an interpolation function f (in 3D) at all vertices of cells. I extract all vertices of cells and then I assign the value to each vertex: if it's in a sphere of radius R, then I assign the value, say 3.91. If it's outside the sphere, then I assign the value 0. I got it running without the error message, but then when I calculated the function value at points that are not vertices, it does not give me either 3.91 or 0. Am I doing something wrong here?

Here is a part of my code:

Extract vertices of all cells, then I export these points and use another software to assign value for each point (say 3.91 for points inside sphere, 0 for points outside)

coor = mesh.coordinates()
numpy.savetxt('meshforE.txt',coor)

I get the values at all vertices and then interpolate this to function f

qvalues2 = numpy.loadtxt('qdata.txt')
V = FunctionSpace(mesh, "CG", 1)
f = Function(V)
f.vector()[:] = qvalues2

then I read points (xp, yp) on z=0 plane and evaluate the funtion at these points

with open('xpdata.txt') as g:
    xp = g.readlines()
print "xp[1]=", xp[1]
with open('ypdata.txt') as h:
    yp = h.readlines()
print "yp[2]=", yp[2]

for i in range(len(xp)):
    g_in[i] = f(xp[i],yp[i],0.0)

Now when I plot gi , it does not look like what it should be, i.e. constant (3.91) inside the circle R=50, and 0 outside. Any help would be appreciated.

Answer 1

You cannot assume that the dofs are numbered in the same way as the vertices. They will be if you add the following parameter:

parameters["reorder_dofs_serial"] = False

but I recommend that you never work with the dofs directly.

Here is simpler and safer way to accomplish what you want to do:

class MyFunction(Expression):
    def eval(self, values, x):
        r2 = x[0]*x[0] + x[1]*x[1]
        if r2 < 0.5*0.5:
            values[0] = 3.91
        else:
            values[1] = 0.0

f = interpolate(MyFunction(), V)

Comments

I want my function f to work with any given data i.e. f cannot be written as a mathematical expression so defining a class "MyFunction" like above won't work. (I set f to be a constant in my question just to check).
parameters["reorder_dofs_serial"] = False
works for me. Thank you.