Using sympy with mayavi

Question

I have a project in which I'd like to plot the results of some SymPy manipulation with mayavi. In particular, I'd like to use some function in contour3d. Here is what I came up with, but it's not working

#!/usr/bin/env python2

from sympy import *
from mayavi.mlab import *
from functools import partial
import numpy as np

def XYZ(f, xp, yp, zp):
  return f.evalf(subs={x:xp,y:yp,z:zp})


x, y, z = symbols('x y z')
r=x**2+y**2+z**2

F=partial(XYZ, r)

X, Y, Z = np.mgrid[-3:3:50j, -3:3:50j, -3:3:50j]

obj=contour3d(X, Y, Z, F)

However, two errors show up (both from sympify)

1.) As written, the error is

ValueError: sequence too large; must be smaller than 32

2.) Decreasing the number of points (to 10 in each direction)

ValueError: negative dimensions are not allowed

Any thoughts on getting something like this to work? I'd think these two packages would go well together.

Answer 1

First of all, I would like to point out that the SymPy team would be very grateful if you contribute a mayavi backend for our plotting module (for the moment we have only a matplotlib backend which is not great for 3d).

The gist of your problem seems to be that numpy and sympy do not cooperate well. Basically, you should not expect (at least for now) that sympy's functions know how to evaluate array inputs (they are not ufuncts as would be said in numpy).

The current solution to this problem is to use lambdify. It has its issues, but it works well enough. A similar, however private function, called experimental_lambdify is used for sympy's own plotting module. See here:

https://github.com/sympy/sympy/blob/master/sympy/plotting/plot.py#L750

Specifically get_meshes(). Or all the other get_something methods. If you are lazy, you can just reuse them by creating a sympy Plot object without actually plotting it in matplotlib:

your_plot_object = plot3d_parametric_surface(expr_x, expr_y, expr_z, range_u, range_v, show=False)
the_single_geometrical_object_in_this_plot = your_plot_object[0]
your_meshes = the_single_geometrical_object_in_this_plot.get_meshes()