NumPy arrays from FEniCS data

Question

I am doing simulations of Navier-Stokes flow on an unstructured grid.

At every timestep I thus have a function u. I know how to export this as an array: u.vector().array() but this is not a very helpful array. In fact, what arrray is it?

What I would like is to somehow get from my mesh and u two arrays:

(1) An array of points on my grid [(x_1,y_1), (x_2, y_2)....]
(2) An array of velocities [(ux_1,uy_1), (ux_1,uy_1)....]

The element at index i in (2) would then be the velocity at the point on the mesh described by the coordinates in place i in (1).

How would I go about creating these arrays from fenics functions?

Thanks in advance.

Here's my attempt, but it is not correct:

arr = u.vector().array()
coor = mesh.coordinates()
vertex_to_dof_map = V.dofmap().vertex_to_dof_map(mesh)

values = list()
for i, dum in enumerate(coor):
    values.append([arr[vertex_to_dof_map[2*i]],arr[vertex_to_dof_map[2*i+1]]])
values = np.array(values)

Answer 1

I found the answer myself.

Here's a full working code:

from dolfin import *
import numpy as np
import pylab as pl

mesh = UnitSquareMesh(2,5)
V = VectorFunctionSpace(mesh, "CG", 1)
u = Function(V)

# Vertex based data
vertex_vector_values = np.zeros(mesh.num_vertices()*2)
vertex_vector_values[::2] = mesh.coordinates().sum(1)
vertex_vector_values[1::2] = 2-mesh.coordinates().sum(1)
vertex_to_dof_map = V.dofmap().vertex_to_dof_map(mesh)

u.vector().set_local(vertex_vector_values[vertex_to_dof_map])

####### NOW GO THE OTHER WAY #####

arr = u.vector().array()
coor = mesh.coordinates()
vtd = V.dofmap().dof_to_vertex_map(mesh)

values = list()
for i, dum in enumerate(coor):
    values.append([arr[vtd[2*i]],arr[vtd[2*i+1]]])
values = np.array(values)

x = list()
y = list()
vx = list()
vy = list()
for i, dum in enumerate(coor):
    print '(%f,%f) -> (%f,%f)' %(coor[i][0], coor[i][1], values[i][0], values[i][1])
    x.append(coor[i][0])
    y.append(coor[i][1])
    vx.append(values[i][0])
    vy.append(values[i][1])

pl.quiver(x,y,vx,vy)
pl.axis([0, 1.3, 0, 1.3])
pl.show()
plot(u, axes=True, interactive=True)

Comments

This works for DOLFIN 1.2. however the vertex_to_dof_map was named wrongly when implemented so it is renamed and moved in the dev version and the upcoming 1.3 version. So instead of:

vertex_to_dof_map = V.dofmap().vertex_to_dof_map(mesh)

you would need to write:

dof_to_vertex_map_values = dof_to_vertex_map(V)

Answer 2

Maybe this code excerpt helps you with your plans. I have once found it in the FEniCS Q&A on launchpad. I assume that V is the function space under consideration.

EDITED: as a follow up of the comments

from dolfin import *
mesh = UnitSquareMesh(3,3)
V = FunctionSpace(mesh, "Lagrange", 1)

for (i, cell) in enumerate(cells(V.mesh())):
        print "Global dofs associated with cell %d: " % i,
        print V.dofmap().cell_dofs(i)
        print "The Dof coordinates:",
        print V.dofmap().tabulate_coordinates(cell)

This will give you the dofs and the corresponding coordinates for each cell. However, since the cells share the dofs, you will have to extract the relevant values.

Comments

This is not super clear... Can we have a complete example starting from a Function to the numpy arrays of coordinates and values, please ?

---

cells(V.mesh) gives an error. Thanks for trying to help though !

---

ok. should be cells(V.mesh())

---

@tlecomte the numpy arrays of the coordinates and values is probably better obtained as in jesper juul's answer.

Answer 3

I offer another answer, for what it is worth.
Define two interpolated values:
x0 = interpolate(Expression("x[0]"),V)
x1 = interpolate(Expression("x[1]"),V)
These two functions have the values of x[0] and x[1] in exactly the same order that your solution, u, has, so that, for i=any legal value, the value u.vector().array()[i] is the value of u at the point (x0.vector().array()[i],x1.vector().array()[i]).