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()

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 $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.


4 Answers 4


Defining a custom Expression and then interpolating that into your function space should do the trick. Check out

from dolfin import *                                                        

mesh = UnitSquareMesh(20, 20)                                               

class CharCircle(Expression):                                               
    def eval(self, value, x):                                               
        xm0 = x[0] - 0.5                                                    
        xm1 = x[1] - 0.5                                                    
        if xm0*xm0 + xm1*xm1 < 0.4**2:                                      
            value[0] = 1.0                                                  
            value[0] = 0.0                                                  

V = FunctionSpace(mesh, 'CG', 1)                                            
u = Function(V)                                                             

  • $\begingroup$ My mesh is 3D though. And I am using external mesh generator and I don't see how your trick would be useful. $\endgroup$
    – Paul S.
    Commented May 15, 2013 at 4:03

I assume that qvalues2 has some computed values at the vertices. You cannot directly assign these to your dof vector as the dofs does not follow vertex numbering.

You could however try:

vertex_to_dof = V.dofmap().vertex_to_dof_map(mesh)
f.vector()[:] = qvalues2[vertex_to_dof]
  • $\begingroup$ Thank you for your answer. Yes, qvalues2 have computed values at vertices. I tried your code but it gave me this error: vertex_to_dof = V.dofmap().vertex_to_dof_map(mesh) AttributeError: 'GenericDofMap' object has no attribute 'vertex_to_dof_map' $\endgroup$
    – Paul S.
    Commented May 9, 2013 at 23:28
  • 1
    $\begingroup$ What version of dolfin do you have? I think that function was added betwee 1.1.0 and 1.2.0. It should have been backported to the 1.1.1 release. $\endgroup$
    – Hake
    Commented May 10, 2013 at 19:12
  • $\begingroup$ I am using Fenics version 1.1.0. Do I need to upgrade and if I do, how do I do that? $\endgroup$
    – Paul S.
    Commented May 14, 2013 at 1:31

You are using piecewise linear functions (V = FunctionSpace(mesh, "CG", 1)), so if you are evaluating the function on a point which is not a vertex, you will get a linear interpolation. You don't say anything about the mesh you are using, but if it's not curvilinear, I would not assume that points on that circle coincide with points on element boundaries (where the interpolant would be constant).

  • $\begingroup$ The mesh I am using is from Gmsh and I converted it. The problem is that for points that are inside the sphere(not on the boundary), my data is all the same (constant=3.91) but when I computed the interpolated function at points that are not vertices, I didn't get 3.91 I got other number between 0-3.91. $\endgroup$
    – Paul S.
    Commented May 9, 2013 at 23:37

If you have the information of your function on all vertices, a representation with Functionspace of continous galerkin degree 1 is only exact at those vertices. Evaluating it at any other point gives the linear interpolation between the points. That's why you observe different values to 3.91 in between values.

If you want that any evaluation within the sphere to be exactly equal to the assigned constant, you need the information of the function on each cell. Maybe by counting a cell to be inside the sphere, when all vertices are inside the sphere. Then you can represent the function with discontinous galerking degree 0. Any evaluation at any point gives than the expected value.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.