A collegue of mine recently tried to visualise strain moduli using matplotlib. If you do not know what that is: Here’s a link if you are interested. https://www.degruyter.com/document/doi/10.1524/zkri.2011.1413/html

It’s a 3D visualisation of the anisotropy of the strain tensor for crystal structures when they are compressed under pressure.

For our purposes it is simply an isosurface in 3D space. You compute values for a 3D meshgrid using a function and then the marching cubes algorithm finds vertex points that map out the isosurface at a given value (let’s say: every vertex point that has an iso value of 0.4).

Additionally, the marching_cubes algorithm finds all faces that belong to the given iso-surface. In the simplest case, there is no anisotropy and the result will be a perfect sphere.

This is the code:

import matplotlib.pyplot as plt
import numpy as np
from skimage import measure


def func(x, y, z):

    Z1111 = 1 
    Z1133 = 1/3 
    Z3333 = 1 

    return (Z1111 * (x**4 + 2 * x**2 * y**2 + y**4) + 6 * Z1133 * (x**2 * z**2 + y**2 * z**2) + Z3333 * z**4 )


a,b,c = -1,1,.1

X, Y, Z = np.mgrid[a:b:c, a:b:c, a:b:c]
vals = func(X,Y,Z)

iso_val = 0.4
verts, faces, murks, _ = measure.marching_cubes(vals, iso_val)

# shift_zero = np.array((a*10,a*10,a*10))
# verts = verts + shift_zero

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot_trisurf(verts[:, 0], verts[:,1], faces, verts[:, 2],
                cmap='jet', lw=1)

ax.plot((0,0),(0,0),(-10,10),'r-',lw=3,zorder=100)
plt.show()

Here’s where the unusual behavior comes in:

When you define the 3D-meshgrid and you use points between -1,1 with an interval of 0.1 for all directions (x,y,z), then the sphere is not centered around (0,0,0), but around (10,10,10). If you uncomment the shift_zero and verts = verts + shift_zero lines, then the offset can be compensated.

Is there something wrong with my script? Does the function for the isosurface include an intrinsic shift? Is this an issue with the implementation of marching_cubes?. I could just compute the center point of the vertices using verts.mean(axis=0) and subtract that from the vertices, but that seems like a cop-out.