Get vorticity from velocity

Hi, I am trying to extract vorticity from the velocity field for fluid simulation. This is my code and extracted field. I am wondering why the vorticity field is so large that it ranges from -160 to 160. Am I doing the computation wrongly?

// We have a known mesh in the file fname
mesh = gen_mesh_from_file(fnames)
nodes = mesh.coordinates()
cell = mesh.cells()

// define function space
V = VectorFunctionSpace(mesh, "Lagrange", 2)
Q = FunctionSpace(mesh, "Lagrange", 1)
u_n = Function(V)

// We have a vector field  u_n in second order VectorFunctionSpace, and we compute vorticity from that 
ux, uy = u_n.split(deepcopy=True)
vortex = ux.dx(1)-uy.dx(0)
vor = Function(Q)
vor = project(vortex,Q)

// We plot the vorticity 
np_solution = vor.compute_vertex_values()
plot_fig = plt.tricontourf(nodes[:,0], nodes[:,1],cell,np_solution,levels=100,

velocity field:


Please follow the guidelines for questions: Read before posting: How do I get my question answered?
and create a Minimal code example that reproduces the issue.

Hi, I have tried to update my question. I can also provide more detailed codes with the mesh file if needed. Thanks a ton.

There are several things that is not clear in your post.

  1. what space is u_n in?
  2. what space is Q?

I have updated the questions!

If u_n is in a second order continuous Lagrange space, then the vorticity is in a Discontinuous Lagrange space of order 1.

Thus your projection to CG-1 is introducing an error between the actual and shown vorticity.

Please use a DG-1 space for the projection, and XDMFFile.write_checkpoint to being able to visualize the DG-1 space in Paraview, as compute_vertex_values only does a CG-1 interpolation of any function.

Hi Dokken, thanks a lot for pointing this out.

I am trying to simulate the fluid dynamics and get the velocity field out of it. Now I am assuming sensors are placed at the node of the mesh. To get the v-field, I am doing this:

V = VectorFunctionSpace(mesh, "Lagrange", 2)
u_1 = Function(V)
// some process to set value for u_1
ux, uy = u_1.split()
u_solution[0, :] = ux.compute_vertex_values()
u_solution[1, :] = uy.compute_vertex_values()  

Do you think it is a correct way of doing this? Here I have u_1 in CG-2 and try to compute_vertex_values() to get the field value.

Another question is that with u_solution, we won’t be able to get u_1 right? As there are many dof are missing. Is it possible to get a u_solution, so that:

V1 = VectorFunctionSpace(mesh, 'P', 1)
u_solution = u_solution.T.reshape(-1)
u1= Function(V1)
d2v_low = dof_to_vertex_map(V1)
u1.vector()[:] = init_u[d2v_low]
V = VectorFunctionSpace(mesh, 'P', 2)
u = project(u1,V)

we have u==u_1.

Why don’t you use the
V.tabulate_dof_coordinates() to get the coordinates of all the degrees of freedom?

As to your second question, you are correct that you cannot recover u_1 by this method, as you have discarded all information about the higher order space when computing the vertex values. You can easily verify this by comparing the underlying arrays of u and u_1 by doing print(u.vector().get_local()-u_1.vector().get_local())

I would still suggest you use XDMFFile.write_checkpoint and Paraview to look at your vorticity.

1 Like

Thanks again Dokken for your help!

It’s a great idea to get V.tabulate_dof_coordinates() to get the coordinates for all the degrees of freedom. Can I also get the value for those nodes? Do you know how should I get this?

One thing I can come up with is to value the field at those coordiantes. Though I am not sure if this is the correct way of doing this.

Hi Dokken, I think I figure this out, I am thinking that I can use


to get the value for each DOF

Yes, the ith entry in u.vector().get_local() corresponds to the ith row of V.tabulate_dof_coordinates()

In particular, see: How to get the solution obatained by fenics as a vector? - #7 by dokken for handling of vector spaces