Can not write vector fields to structured grid

[karpathyan]

Hi, thanks for this nice tool. But When I try to use this to write vector files vector fields to the VTS file, I'm getting an error. Can you please look into this? The code I used and the error obtained are given below.

code (this is a modified version of the structured.py example file)

import os
from VTKwrite.interface import structuredToVTK
import numpy as np
import random as rnd

FILE_PATH = "./structured"
def clean():
    try:
        os.remove(FILE_PATH + ".vts")
    except:
        pass

def run():
    print("Running structured...")

    # Dimensions
    nx, ny, nz = 6, 6, 2
    lx, ly, lz = 1.0, 1.0, 1.0
    dx, dy, dz = lx/nx, ly/ny, lz/nz

    ncells = nx * ny * nz
    npoints = (nx + 1) * (ny + 1) * (nz + 1)

    # Coordinates
    X = np.arange(0, lx + 0.1*dx, dx, dtype='float64')
    Y = np.arange(0, ly + 0.1*dy, dy, dtype='float64')
    Z = np.arange(0, lz + 0.1*dz, dz, dtype='float64')

    x = np.zeros((nx + 1, ny + 1, nz + 1))
    y = np.zeros((nx + 1, ny + 1, nz + 1))
    z = np.zeros((nx + 1, ny + 1, nz + 1))

    # Variables
    pressure = np.random.rand(ncells).reshape( (nx, ny, nz))
    temp = np.random.rand(npoints).reshape( (nx + 1, ny + 1, nz + 1))
    my_vec = np.random.rand(nx, ny, nz, 3)
    # data defined on cells
    all_cell_data = []
    print (np.shape(np.array(my_vec).reshape(nx,ny,nz,3)))
    print ((nx, ny, nz))
    all_cell_data.append(["my_vec","vectors", my_vec ])

    # data defined at points
    all_point_data = []
    all_point_data.append(["temp", "scalars", temp])

    comments = [ "comment 1", "comment 2" ]
    structuredToVTK(FILE_PATH, x, y, z, all_cell_data = all_cell_data, all_point_data = all_point_data, comments = comments)

if __name__ == "__main__":
    run()

Error obtained

Traceback (most recent call last):
  File "/PyVTK_write/VTKwrite-main/examples/Structured/structured.py", line 59, in <module>
    run()
  File "PyVTK_write/VTKwrite-main/examples/Structured/structured.py", line 56, in run
    structuredToVTK(FILE_PATH, x, y, z, all_cell_data = all_cell_data, all_point_data = all_point_data, comments = comments)
  File "/home/.local/lib/python3.10/site-packages/VTKwrite-0.0.2-py3.10.egg/VTKwrite/interface.py", line 304, in structuredToVTK
AssertionError

[walkersw]

Thank you for this information. I have now fixed the code. Re-download everything, re-install, and look at the new example file: structured_ex1.py

[karpathyan]

Hi, Thanks for the quick fix the new example is working fine. However, I find the implementation a little confusing. Can you please explain the working to me? For example, in structured_ex1.py, we have a 6x6x2 grid. So the vector field should have a shape of (6,6,2,3) right? However, in line 93 when we are providing the Vec_cell to the append function it has the shape (216,). So basically, you are flattening the (6,6,2,3) Ndimessional array before appending it to the grid. So how can we keep track of which vectors are going into which cell?

So for context, I'm planning to use your code for creating VTK files from a series of Finite difference vector data, That is I have my vector field in the following format,

px py pz   Vx Vy Vz

Where px, py, pz are the cartesian coordinates of the cells and Vx, Vy, Vz are the vector components of the cells.

So if I have a 10x10x2 grid, I will have 200 cells and 200 vectors. I have the cartesian coordinates of the cells and their vector information. How can I implement this using your code to create a vts file?

I Thank you for your time.

[walkersw]

You should be able to generalize the example. For example, change nx and ny to 10 in your case above. The example code even samples a given scalar and vector field.

[karpathyan]

Can you please comment on this?

For example, in structured_ex1.py, we have a 6x6x2 grid. So the vector field should have a shape of (6,6,2,3) right? However, in line 93 when we are providing the Vec_cell to the append function it has the shape (216,). So basically, you are flattening the (6,6,2,3) Ndimessional array before appending it to the grid. So how can we keep track of which vectors are going into which cell?

[walkersw]

Please read the example code. For example, Vec_cell is initialized with

Vec_cell = np.zeros((3, nx, ny, nz))

which means it has shape (3, 6, 6, 2), because nx=6, ny=6, and nz=2. you then fill that array (see the example code), and then you flatten:

Vec_cell = Vec_cell.reshape(3*ncells, order='F')

where ncells = nxnynz.

[karpathyan]

Thanks for your help, I think I was not able to express my doubt clearly. My question is, for example if we have a 3x3x1 grid: the first cell will be [0,0,0], What will be the second cell? is it [1,0,0] or [0,0,1] or [0,1,0] ? How are the cells in your array arranged? This is important because we have to match the vector array with this right?

[walkersw]

You will need to do some testing on your part to check. It sounds like you want to implement the grid slightly different from how I have. If you do it the way I have, then you don't need to know the details. Otherwise, you will need to run some examples on your end to figure it out.

You could try flattening (in the same way) the x, y, z grid points, or xm, ym, zm, to check the ordering. Either way you are going to have to test, i.e. do small examples and print out the coordinates, etc.