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=", xp with open('ypdata.txt') as h: yp = h.readlines() print "yp=", yp for i in range(len(xp)): g_in[i] = f(xp[i],yp[i],0.0)
Now when I plot $g_i$ , 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.