What is the best way to evaluate a Function in a list of coordinates?

Question

I am trying to evaluate a Function coming from a solve call in FEniCS 1.5.0 and tried this in the first place:

NUM_ELEM = 10
mesh = UnitSquareMesh(NUM_ELEM, NUM_ELEM)
# [...]
vals = np.empty((NUM_ELEM + 1) * (NUM_ELEM + 1), dtype=float)
u.eval(vals, mesh.coordinates())

but got a TypeError: contiguous numpy array of 'double' expected. Make sure that the numpy array is contiguous, with 1 dimension, and uses dtype=float_. I have tried also this combination but fails with the same exception:

for v in vertices(mesh):
    u.eval(vals[:, None][v.index()], mesh.coordinates()[v.index()])

What is the best way to evaluate a Function in an array? I would like a solution working in a general case, i.e. points not coming from the mesh, but from np.meshgrid for example. I am aware I could always write a for loop but I am looking for a straightforward approach.

Comments

You are using python type float, but you should use type numpy.float_.

---

I tried several combinations but still doesn't work:

In [7]: vals = np.empty(((NUM_ELEM + 1) * (NUM_ELEM + 1), 2), dtype=np.float_)

In [8]: vals.shape
Out[8]: (121, 2)

In [9]: vals.dtype
Out[9]: dtype('float64')

In [10]: u.eval(vals, mesh.coordinates())
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-10-5b538805dc98> in <module>()
----> 1 u.eval(vals, mesh.coordinates())

TypeError: contiguous numpy array of 'double' expected. Make sure that the numpy array is contiguous, with 1 dimension, and uses dtype=float_.

In [11]: vals = np.empty(((NUM_ELEM + 1) * (NUM_ELEM + 1),), dtype=np.float_)

In [12]: vals.shape
Out[12]: (121,)

In [13]: vals.dtype
Out[13]: dtype('float64')

In [14]: u.eval(vals, mesh.coordinates())
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-14-5b538805dc98> in <module>()
----> 1 u.eval(vals, mesh.coordinates())

TypeError: contiguous numpy array of 'double' expected. Make sure that the numpy array is contiguous, with 1 dimension, and uses dtype=float_.

In [15]: mesh.coordinates().shape
Out[15]: (121, 2)

Answer 1

There are solutions to do just this in fenicstools. See, e.g., probes

Different solution:

for v in vertices(mesh):
    vals[v.index()] = u_(*mesh.coordinates()[v.index()])

Comments

This is very useful, thanks! If there is no other way I will accept this answer.

---

There are always different ways:-) Updated the answer.

---

My hope was in avoiding the for loop :) Accepting the answer!

---

You can interpolate on a RectangleMesh.

---

Nice idea, perhaps I could transform a meshgrid into a RectangleMesh.

---

I am pretty sure the ordering of a CG1 field on a RectangleMesh is identical to what you would get for a MATLAB style meshgrid.

Answer 2

I believe the error message is not related to vals, but to mesh.coordinates(), which is also supposed to 1-dimension.
This should work,

u.eval(vals, mesh.coordinates().flatten())