Search code examples
pythoncoordinate-systemspyvista

How to convert a cartesian problem in a cylindrical problem?

I display a gyroid structure (TPMS) in a cartesian system using Pyvista. I try now to display the structure in cylindrical coordinates. Pyvista displays something cylindrical indeed but it seems that the unit cell length is not uniform (while there is no reason to change this my parameter "a" being steady. This change seems to appear especially along z but I don't understand why (see image).

I have this: enter image description here Here is a part of my code.

Thank you for your help.

import pyvista as pv
import numpy as np
from numpy import cos, sin, pi
from random import uniform

lattice_par = 1.0  # Unit cell length
a = (2*pi)/lattice_par

res = 200j
r, theta, z = np.mgrid[0:2:res, 0:2*pi:res, 0:4:res]
# consider using non-equidistant r for uniformity

def GyroidCyl(r, theta, z, b=0.8):
    return (sin(a*(r*cos(theta) - 1))*cos(a*(r*sin(theta) - 1))
            + sin(a*(r*sin(theta) - 1))*cos(a*(z - 1))
            + sin(a*(z - 1))*cos(a*(r*cos(theta) - 1))
            - b)

vol3 = GyroidCyl(r, theta, z)

# compute Cartesian coordinates for grid points
x = r * cos(theta)
y = r * sin(theta)

grid = pv.StructuredGrid(x, y, z)
grid["vol3"] = vol3.flatten()
contours3 = grid.contour([0])  # Isosurface = 0

pv.set_plot_theme('document')
p = pv.Plotter()
p.add_mesh(contours3, scalars=contours3.points[:, 2], show_scalar_bar=False, interpolate_before_map=True,
           show_edges=False, smooth_shading=False, render=True)
p.show_axes_all()
p.add_floor()

p.show_grid()
p.add_title('Gyroid in cylindrical coordinates')
p.add_text('Volume Fraction Parameter = ' + str(b))
p.show(window_size=[2040, 1500])
Solution

  • So you've noted in comments that you're trying to replicate something like the strategy explained in this paper. What they do is take a regular gyroid unit cell, and then transform it to build a cylindrical shell. If igloos were cylindrical, then a gyroid cell would be a single piece of snow brick. Put them next to one another and stack them in a column, and you get a cylinder.

    Since I can't use figures from the paper we'll have to recreate some ourselves. So you have to start from a regular gyroid defined by the implicit function

    cos(x) sin(y) + cos(y) sin(z) + cos(z) sin(x) = 0

    (or some variation thereof). Here's how a single unit cell looks:

    import pyvista as pv
import numpy as np

res = 100j
a = 2*np.pi
x, y, z = np.mgrid[0:a:res, 0:a:res, 0:a:res]

def Gyroid(x, y, z):
    return np.cos(x)*np.sin(y) + np.cos(y)*np.sin(z) + np.cos(z)*np.sin(x)

# compute implicit function
fun_values = Gyroid(x, y, z)

# create grid for contouring
grid = pv.StructuredGrid(x, y, z)
grid["vol3"] = fun_values.ravel('F')
contours3 = grid.contour([0])  # isosurface for 0

# plot the contour, i.e. the gyroid
pv.set_plot_theme('document')
plotter = pv.Plotter()
plotter.add_mesh(contours3, scalars=contours3.points[:, -1],
                 show_scalar_bar=False)
plotter.add_bounding_box()
plotter.enable_terrain_style()
plotter.show_axes()
plotter.show()

    single unit cell of the conventional gyroid lattice

    Using the "unit cell" term implies there's an underlying infinite lattice, which can be built by stacking these (rectangular) unit cells neatly next to one another. With some imagination we can convince ourselves that this is true. Or we can look at the formula and note that due to the trigonometric functions the function is periodic in x, y and z, with period 2*pi. This also tells us that we can transform the unit cell to have arbitrary rectangular dimensions by introducing lattice parameters a, b and c:

    cos(kx x) sin(ky y) + cos(ky y) sin(kz z) + cos(kz z) sin(kx x) = 0, where
kx = 2 pi/a
ky = 2 pi/b
kz = 2 pi/c

    (These kx, ky and kz quantities are called wave vectors in solid state physics.)

    The corresponding change only affects the header:

    res = 100j
a, b, c = lattice_params = 1, 2, 3
kx, ky, kz = [2*np.pi/lattice_param for lattice_param in lattice_params]
x, y, z = np.mgrid[0:a:res, 0:b:res, 0:c:res]

def Gyroid(x, y, z):
    return (  np.cos(kx*x)*np.sin(ky*y)
            + np.cos(ky*y)*np.sin(kz*z)
            + np.cos(kz*z)*np.sin(kx*x))

    gyroid with (1, 3, 5) lattice parameters

    This is where we start. What we have to do is take this unit cell, bend it so that it corresponds to a 30-degree circular arc on a cylinder, and stack the cylinder using this unit. According to the paper, they used 12 unit cells to create a circle in a plane (hence the 30-degree magic number), and stacked three such circular bands on top of each other to build the cylinder.

    The actual mapping is also fairly clearly explained in the paper. Whereas your original x, y and z parameters of the function essentially interpolated between [0, a], [0, b] and [0, c], respectively, in the new setup x interpolates in the radius range [r1, r2], y interpolates in the angular range [0, pi/6] and z is just z. (In the paper x and y seem to be reversed with respect to this convention, but this shouldn't matter. If it matters, that's left as an exercise to the reader.)

    So what we need to do is more or less keep the current grid points, but transform the corresponding x, y and z grid points so that they lie on a cylinder instead. Here's one take:

    import pyvista as pv
import numpy as np

res = 100j
a, b, c = lattice_params = 1, 1, 1
kx, ky, kz = [2*np.pi/lattice_param for lattice_param in lattice_params]
r_aux, phi, z = np.mgrid[0:a:res, 0:b:res, 0:3*c:res]

# convert r_aux range to actual radii
r1, r2 = 1.5, 2
r = r2/a*r_aux + r1/a*(1 - r_aux)

def Gyroid(x, y, z):
    return (  np.cos(kx*x)*np.sin(ky*y)
            + np.cos(ky*y)*np.sin(kz*z)
            + np.cos(kz*z)*np.sin(kx*x))

# compute data for cylindrical gyroid
# r_aux is x, phi / 12 is y and z is z
fun_values = Gyroid(r_aux, phi * 12, z)

# compute Cartesian coordinates for grid points
x = r * np.cos(phi*ky)
y = r * np.sin(phi*ky)
grid = pv.StructuredGrid(x, y, z)
grid["vol3"] = fun_values.ravel('F')
contours3 = grid.contour([0])

# plot cylindrical gyroid
pv.set_plot_theme('document')
plotter = pv.Plotter()
plotter.add_mesh(contours3, scalars=contours3.points[:, -1],
                 show_scalar_bar=False)
plotter.add_bounding_box()
plotter.show_axes()
plotter.enable_terrain_style()
plotter.show()

    cylindrical gyroid

    If you want to look at a single transformed unit cell in the cylindrical setting, use a single domain of phi and z for the function and only convert to 1/12 a full circle for the grid points:

    fun_values = Gyroid(r_aux, phi, z/3)

# compute Cartesian coordinates for grid points
x = r * np.cos(phi*ky/12)
y = r * np.sin(phi*ky/12)
grid = pv.StructuredGrid(x, y, z/3)

    But it's not easy to see the curvature in the (no longer a) unit cell:

    transformed unit cell