 # 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,cmap=cmocean.cm.thermal)

``````

velocity field:

vorticity: 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

``````u.vector().get_local()
``````

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