Hex mesh from finite difference voxels

[tloudon]

Hi all,

I’m working on generating a hexahedral mesh in Cubit from voxelized data produced by a Finite Difference Time Domain (FDTD) simulation.

Each voxel is axis-aligned and defined by its center coordinates and dimensions. The goal is to create one HEX8 element per voxel, ideally resulting in a merged and compact mesh that can be exported to Exodus format.

I’ve written a minimal working example (MWE) that reads in the voxel data and:

• Creates one brick per voxel
• Moves each brick to the correct location
• Meshes each brick as a single hexahedral element

This works for small cases, but performance is a major concern. For example:

• Meshing 70,000 voxels takes approximately 30 minutes
• My production use cases involve between 10 million and 100 million voxels

To avoid further slowdown, I’m currently skipping the geometry merge step, though I would prefer to merge the result to reduce file size. I haven’t yet used Cubit’s parallel features, but this seems like a good use case for them. I also have significant computational resources to bring to bear for this task.

Any suggestions or ideas to improve performance are very welcome.

#Import packages
import numpy as np
import cubit
from timeit import time

#Set file names
voxel_data_file = 'voxel_data.txt.gz'
exo_file = 'voxel_mesh.exo'

#Initialize cubit
cubit.init(['cubit', '-nojournal', '-nographics'])
cubit.reset()

#Some default settings for improved performance
cubit.cmd('echo off')
cubit.cmd('undo off')
cubit.cmd('warning off')
cubit.cmd('info off')
cubit.cmd('set default autosize off')

#Load voxel centers
xcm, ycm, zcm, Lx, Ly, Lz = np.loadtxt(voxel_data_file)
nvoxel = int(xcm.shape[0])
vlength = len(str(nvoxel))

#Create bricks for each voxel
t0 = time.time()
for k in range(nvoxel):
    cubit.cmd(f'brick x {Lx[k]} y {Ly[k]} z {Lz[k]}')
    voxel_id = cubit.get_last_id('volume')
    cubit.cmd(f'volume {voxel_id} move x {xcm[k]} y {ycm[k]} z {zcm[k]}')
    cubit.cmd(f'block 1 vol {voxel_id}')
    cubit.cmd(f'vol {voxel_id} interval 1')
    cubit.cmd(f'vol {voxel_id} scheme map')
    cubit.cmd(f'mesh vol {voxel_id}')
    print(f'Voxel {k+1:0{vlength}d} of {nvoxel} created at time {time.time()-t0:.2f}', flush=True)
assert cubit.get_element_count() == nvoxel, "Something is off with the mesh."

#Add to block
cubit.cmd('block 1 element HEX8')

#Save
cubit.cmd(f'export genesis "{exo_file}" overwrite')
print(f'\nVoxel mesh created in {time.time()-t0:.2f} seconds.', flush=True)

voxel_data.txt.gz (49.9 KB)

[Norbert_Hofbauer]

Hi @tloudon,
should all your voxels have the same dimension?

When i look at Lx, Ly, Lz i would say a constant 4 would be good unless you say the dimensions will vary depending on the voxel data.

[tloudon]

In general the voxels will be nonuniform even if they are (mostly) uniform for this simple example problem. Thanks.

[tloudon]

To anyone interested, here is what chatGPT and I came up with. The data is in a slightly different format than I posted previously. The coords array is the 3D coordinates of the 8 corners of each voxel. This way able to handle a FDTD mesh with 62.8M voxels in just over an hour.

#Import packages
import numpy as np
import meshio
from scipy.spatial import cKDTree
from netCDF4 import Dataset

# Parameters
voxel_data_file='voxel_data.txt.gz'
outfile = 'voxel_mesh.exo'
tolerance = 1e-1 #Consider points within a distance < tolerance equal

# Load voxel data
data = np.loadtxt(voxel_data_file, skiprows=1, delimiter=',')
coords = data[:, 2:5] #(x, y, z) coordinates of each voxel corner
og_rows = len(coords)

# Build KDTree
tree = cKDTree(coords)

# Get connected components using a union-find (disjoint-set) approach
n = len(coords)
parent = np.arange(n)
def find(i):
    while parent[i] != i:
        parent[i] = parent[parent[i]]
        i = parent[i]
    return i

# Query all pairs within tolerance and group them
pairs = tree.query_pairs(r=tolerance)
for i, j in pairs:
    pi, pj = find(i), find(j)
    if pi != pj:
        parent[pj] = pi

# Assign group ID to each point
group_ids = np.array([find(i) for i in range(n)])

# group_ids: shape (N,), from union-find
max_gid = group_ids.max()
relabel = np.full(max_gid + 1, -1, dtype=np.int32)
used_gids = np.flatnonzero(np.bincount(group_ids))  # Unique group IDs (very fast)
relabel[used_gids] = np.arange(len(used_gids))
inverse_indices = relabel[group_ids]  # Shape (N,), values in [0, num_unique)
num_unique = inverse_indices.max() + 1

# Output
print(f"Original coordinates: {og_rows}", flush=True)
print(f"Unique coordinates:   {num_unique}", flush=True)
print(f"Duplicates removed:   {og_rows - num_unique}", flush=True)

# Compute canonical coords using fast binning
canonical_coords = np.zeros((num_unique, 3))
for i in range(3):
    canonical_coords[:, i] = np.bincount(inverse_indices, weights=coords[:, i], minlength=num_unique)
counts = np.bincount(inverse_indices, minlength=num_unique)
canonical_coords /= counts[:, None]

# Final connectivity
connectivity = inverse_indices.reshape((-1, 8)).astype(np.int32)

#Create meshio mesh
cells = [("hexahedron", connectivity)]
meshio.write_points_cells(outfile, canonical_coords, cells)

#Shift block 0 -> 1
with Dataset(outfile, 'r+') as ds:
    ebids = ds.variables['eb_prop1'][:]
    ds.variables['eb_prop1'][:] = 1

[Norbert_Hofbauer]

Hi @tloudon,
i now tried a few approaches which generally would work, but are way to slow to handle millions of voxels which are not of the same size.

Just for example. Creating a high number of single hexes is quite slow using cubit commands, so i would rather create a abaqus file that can be imported. With equivalence node all tolerance 1e-1 i can merge all nodes within a tolerance. But unfortunately this is quite slow for a big number of nodes and it gets slower the higher the node count is. Even when i try to split the number of nodes there is not really much speed up.

Your solution is definitely working faster, so you should stick to it. I will report the performance issue using equivalence node all.

#!python

import sys

# Append path to Cubit module to sys.path
sys.path.append('/opt/Coreform-Cubit-2025.3/bin')

import numpy
import cubit

cubit.init(['cubit', '-nojournal', '-nographics'])
cubit.reset()

def write_nodes(f,k,x, y, z, dx, dy, dz):
 f.write(f"{1+8*k}, {x+dx/2}, {y-dy/2}, {z+dz/2}\n")   
 f.write(f"{2+8*k}, {x+dx/2}, {y-dy/2}, {z-dz/2}\n")   
 f.write(f"{3+8*k}, {x-dx/2}, {y-dy/2}, {z-dz/2}\n")   
 f.write(f"{4+8*k}, {x-dx/2}, {y-dy/2}, {z+dz/2}\n")   
 f.write(f"{5+8*k}, {x+dx/2}, {y+dy/2}, {z+dz/2}\n")   
 f.write(f"{6+8*k}, {x+dx/2}, {y+dy/2}, {z-dz/2}\n")   
 f.write(f"{7+8*k}, {x-dx/2}, {y+dy/2}, {z-dz/2}\n")   
 f.write(f"{8+8*k}, {x-dx/2}, {y+dy/2}, {z+dz/2}\n")   

def write_hex(f,k):
 f.write(f"{1+k}, {1+8*k}, {2+8*k}, {3+8*k}, {4+8*k}, {5+8*k}, {6+8*k}, {7+8*k}, {8+8*k}\n")   

#Set file names
voxel_data_file = 'voxel_data.txt.gz'
abq_inp = 'test.inp'
exo_file = 'voxel_mesh.exo'

cubit.silent_cmd('echo off')
cubit.silent_cmd('undo off')
cubit.silent_cmd('warning off')
cubit.silent_cmd('info off')

#Load voxel centers
xcm, ycm, zcm, Lx, Ly, Lz = numpy.loadtxt(voxel_data_file)

#nodes
nvoxel = int(xcm.shape[0])
i = 0

with open(abq_inp, "w") as f:
 f.write("*NODE, NSET=ALLNODES\n")
 for k in range(nvoxel):
  write_nodes(f,k,xcm[k], ycm[k], zcm[k], Lx[k], Ly[k], Lz[k])
  i = i +1
  if i == 1000:
   status=cubit.silent_cmd('echo on')
   print(f"NODES: processed {k+1} from {nvoxel} voxels")
   status=cubit.silent_cmd('echo off')
   i=0
  if k==5000:
   break

#hex connectivity
with open(abq_inp, "a") as f:
 f.write("*ELEMENT, TYPE=C3D8, ELSET=EB1\n")
 for k in range(nvoxel):
  write_hex(f,k)
  i = i +1
  if i == 1000:
   status=cubit.silent_cmd('echo on')
   print(f"HEX: processed {k+1} from {nvoxel} voxels")
   status=cubit.silent_cmd('echo off')
   i=0
  if k==5000:
   break

cubit.silent_cmd('echo on')

cubit.silent_cmd('import abaqus  "/home/user/coreform/test.inp"')
#equivalence node all tolerance 1e-1 preview
#equivalence node all tolerance 1e-1

#!python
xmin = -100
xmax = 100
dx = 4
print(f"xmin = {xmin},xmax = {xmax}")
for x in range(xmin,xmax,dx):
 ids = cubit.parse_cubit_list("node",f"all with x_coord>={x-dx/2} and x_coord<={x+dx/2}")
 if len(ids)>0:
  status = cubit.silent_cmd(f"equivalence node {' '.join(str(id) for id in ids)} tolerance 1e-1")
  print(f"x = {x}, equivalence {len(ids)} nodes")

#cubit.cmd(f'export genesis "{exo_file}" overwrite')

#!cubit
cubit.silent_cmd('echo on')
cubit.silent_cmd('undo on')
cubit.silent_cmd('warning on')
cubit.silent_cmd('info on')