This Jupyter Notebook explores voxel-based simulations for physical systems. It discusses optimizing voxel tessellation algorithms, including the considerations for using cubic or truncated octahedral voxel shapes. It also entails practical implementation of these concepts in Python, utilizing numpy and matplotlib for operations like defining 3D kernels for neighbor interactions in voxel grids, computing temperature differences for a modeled sphere, and visualizing these voxel grids. At its core, the notebook serves as a tutorial or example code for conducting physics simulations or studying material properties at the voxel level, including aspects like energy conservation, pressure, and temperature variations within the system.

- Should I consider changing the tesselation from cubic to truncated Octahedral - This would make some things easier (like all neighbors would share area with all other neighbors). However, some things would be harder, for example the area of each face is not constant so I have to keep track then of what kind of face makes contact. Further, I think the seperation between grid points becomes a bit harder to work with)

snippet.python

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.colors import Normalize
from matplotlib.cm import ScalarMappable
from scipy.ndimage import convolve, correlate
from mplEasyAnimate import animation
from matplotlib.cm import viridis
from tqdm.notebook import tqdm
plt.style.use(pubStyle)

snippet.python

def make_3d_kernels():
    # Central kernel
    ck = np.zeros((3, 3, 3))
    ck[1, 1, 1] = 1
 
    # Initialize kernels array to hold all 26 kernels
    kernels = np.zeros((26, 3, 3, 3))
 
    # Counter for filling the kernels array
    counter = 0
 
    # Generate face neighbor kernels
    for axis in range(3):  # Loop through x, y, z axes
        for delta in [-1, 1]:  # Negative and positive directions
            k = np.copy(ck)
            k[1 + delta*(axis == 0), 1 + delta*(axis == 1), 1 + delta*(axis == 2)] = -1
            kernels[counter, :, :, :] = k
            counter += 1
 
    # Generate edge neighbor kernels
    for x in [-1, 1]:
        for y in [-1, 1]:
            k = np.copy(ck)
            k[1 + x, 1 + y, 0] = -1/2
            kernels[counter, :, :, :] = k
            counter += 1
 
            k = np.copy(ck)
            k[1 + x, 0, 1 + y] = -1/2
            kernels[counter, :, :, :] = k
            counter += 1
 
            k = np.copy(ck)
            k[0, 1 + x, 1 + y] = -1/2
            kernels[counter, :, :, :] = k
            counter += 1
 
    # Generate corner neighbor kernels
    for x in [-1, 1]:
        for y in [-1, 1]:
            for z in [-1, 1]:
                k = np.copy(ck)
                k[1 + x, 1 + y, 1 + z] = -1/3
                kernels[counter, :, :, :] = k
                counter += 1
 
    return kernels

snippet.python

def pretty_print_3d_array(array3D):
    # Determine the maximum number width for formatting
    maxNumWidth = max(len(str(abs(x))) for x in array3D.flatten() if x != 0)
 
    # Iterate through each layer
    for layerIndex in range(array3D.shape[0]):
        for rowIndex in range(array3D.shape[1]):
            for colIndex in range(array3D.shape[2]):
                element = array3D[layerIndex, rowIndex, colIndex]
                # Using format specifiers to align the output
                print(f"{element:{maxNumWidth}f}" if element != 0 else " " * maxNumWidth, end=" ")
            # Move to the next row
            print()
        # Add an empty line between layers for better readability
        print()

snippet.python

class EnergyConservationError(Exception):
    def __init__(self, msg):
        self.message = msg

snippet.python

class TimeEvolver:
    def __init__(self, sphere, times):
        self.sphere = sphere
        self._tracker = np.zeros_like(shape=(*self.sphere.temperature.shape, 4, len(times)))
        self.times = times
        self.dts = np.gradient(self.times)
 
    def evolve(self):
        for tID, dt in enumerate(self.dts):
            self.sphere.timestep(dt)

snippet.python

class CallbackNDArray(np.ndarray):
    def __new__(cls, input_array, callback=None, *args, **kwargs):
        obj = np.asarray(input_array).view(cls)
        obj.callback = callback
        return obj
 
    def __setitem__(self, key, value):
        super().__setitem__(key, value)
        if hasattr(self, 'callback') and self.callback:
            self.callback()

snippet.python

class VoxelSphere:
    mH = 1.00784  
    mHe = 4.002602  
    R = 8.3145e7
    R = 1
 
    _3DDiffKernels = make_3d_kernels()
 
    def __init__(self, radius, resolution, t0=0, X=0.75, Y=0.25, Z=0, D=1, k=1):
        self._radius = radius
        self._resolution = resolution
        self._t = t0
        self._X, self._Y, self._Z = X, Y, Z
        self._effectiveMolarMass = self._X * self.mH + self._Y * self.mHe
 
        self._create_voxel_sphere()
 
        self._D = D
        self._k = k
 
    @property
    def radius(self):
        return self._radius
 
    @property
    def resolution(self):
        return self._resolution
 
    def _create_voxel_sphere(self):
        gridSize = np.ceil(self.radius * self.resolution) * 2 + 1
 
        self.x = np.linspace(-self.radius, self.radius, int(gridSize))
        self.y = np.linspace(-self.radius, self.radius, int(gridSize))
        self.z = np.linspace(-self.radius, self.radius, int(gridSize))
        xGrid, yGrid, zGrid = np.meshgrid(self.x, self.y, self.z, indexing='ij')
        self._dr = np.linalg.norm(np.array([xGrid[0,0,0] - xGrid[0,1,0], yGrid[0,0,0] - yGrid[1,0,0], zGrid[0,0,0] - zGrid[0,0,1]]))
        self._voxelVolume = np.abs((self.x[1] - self.x[0]) * (self.y[1] - self.y[0]) * (self.z[1] - self.z[0]))
        self._voxelFaceArea = np.abs((self.x[1] - self.x[0]) * (self.y[1] - self.y[0]))
 
        self._temperature = CallbackNDArray(np.ones_like(xGrid), callback=self._recompute_energy)
        self._density = CallbackNDArray(np.ones_like(xGrid), callback=self._recompute_energy)
        self._pressure = CallbackNDArray(np.ones_like(xGrid), callback=self._recompute_energy)
 
        distanceFromCenter = np.sqrt(xGrid**2 + yGrid**2 + zGrid**2)
        self._voxelSphere = distanceFromCenter <= self.radius
        self._recompute_energy()        
 
    def _recompute_energy(self):
        self._n = (self._pressure * self._voxelVolume)/self._effectiveMolarMass
        self._energy = (3/2) * self._n * self.R * self._temperature
        self._energy[~self._voxelSphere] = 0
 
    def _recompute_temperature(self):
        self._n = (self._pressure * self._voxelVolume)/self._effectiveMolarMass
        self._temperature = (2/3) * self._energy/(self._n*self.R)
        self._temperature[~self._voxelSphere] = 0
 
    @property
    def temperature(self):
        return self._temperature
 
    def _temperature_differences(self):
        diffs = np.zeros(shape=(26, *self._temperature.shape))
        for kernelID, kernel in enumerate(self._3DDiffKernels):
            correlation = correlate(self._temperature, kernel, mode='constant', cval=0.0)
            diffs[kernelID] = correlation
        return diffs
 
    def timestep(self, dt, ect = 0.01):
        totalEnergyInit = self._energy.sum()
        dTs = self._temperature_differences()
        Q = np.zeros_like(dTs)
        for orientationID, dT in enumerate(dTs): 
            orientationQ = -self._k*self._voxelFaceArea*dT/(self._dr)
        self._energy = self._energy + dE
        pretty_print_3d_array(self._energy)
        totalEnergyFinal = self._energy.sum()
        fchange = abs((totalEnergyFinal - totalEnergyInit)/totalEnergyInit)
        if fchange > ect:
            raise EnergyConservationError(f"Total Energy varied more than allowed minimum fraction ({ect*100:0.3f}%) {fchange*100:0.3f}%")
        self._recompute_temperature()
 
    def plot_voxel_region(self, voxelSphere, voxelColors='red', **kwargs):
        fig, ax = plt.subplots(1, 1, figsize=(10, 7), subplot_kw=dict(projection='3d'))
 
        filled = np.zeros(voxelSphere.shape, dtype=bool)
        filled[np.where(voxelSphere)] = True
        ax.voxels(filled, facecolors=voxelColors, edgecolor='k', **kwargs)
 
        ax.set_xlabel('X Axis')
        ax.set_ylabel('Y Axis')
        ax.set_zlabel('Z Axis')
 
        return fig, ax
 
    def _mk_voxel_color_array(self, array):
        Norm = Normalize(vmin=array.min(), vmax=array.max())
        colors = viridis(Norm(array).flatten()).reshape(*array.shape, 4)
        return colors
 
    def _pcparam_to_color(self, pcparam):
        if pcparam == None:
            colors = np.zeros(shape=(*self._voxelSphere.shape, 4))
            colors[:,:,:] = [1,1,1,1]
        elif pcparam == 'temperature':
            colors = self._mk_voxel_color_array(self._temperature)
        elif pcparam == 'pressure':
            colors = self._mk_voxel_color_array(self._pressure)
        elif pcparam == 'energy':
            colors = self._mk_voxel_color_array(self._energy)
        return colors
 
    def plot_voxel_sphere(self, pcparam=None, cmap='viridis', **kwargs):
        colors = self._pcparam_to_color(pcparam)
        return self.plot_voxel_region(self._voxelSphere, voxelColors=colors, **kwargs)

snippet.python

snippet.python

sphere = VoxelSphere(1, 1)
sphere.timestep(1)

                                                         
                             0.853956                    
                                                         

                             0.853956                    
          0.853956           0.853956           0.853956 
                             0.853956                    

                                                         
                             0.853956                    
                                                         

================================
                                                      
                           4.269781                   
                                                      

                           4.269781                   
         4.269781          5.123737          4.269781 
                           4.269781                   

                                                      
                           4.269781                   
                                                      

================================
                                                      
                           5.123737                   
                                                      

                           5.123737                   
         5.123737          5.977693          5.123737 
                           5.123737                   

                                                      
                           5.123737                   
                                                      




---------------------------------------------------------------------------

EnergyConservationError                   Traceback (most recent call last)

Cell In[236], line 2
      1 sphere = VoxelSphere(1, 1)
----> 2 sphere.timestep(1)


Cell In[235], line 74, in VoxelSphere.timestep(self, dt, ect)
     72 fchange = abs((totalEnergyFinal - totalEnergyInit)/totalEnergyInit)
     73 if fchange > ect:
---> 74     raise EnergyConservationError(f"Total Energy varied more than allowed minimum fraction ({ect*100:0.3f}%) {fchange*100:0.3f}%")
     75 self._recompute_temperature()


EnergyConservationError: Total Energy varied more than allowed minimum fraction (1.000%) 514.286%

snippet.python

a = np.zeros((10, 10))
a[5, :] = 5
a[6, :] = 6
k = np.zeros((3,3))
k[1,1] = 1
k[0,1] = -1
a

array([[0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [5., 5., 5., 5., 5., 5., 5., 5., 5., 5.],
       [6., 6., 6., 6., 6., 6., 6., 6., 6., 6.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.]])

snippet.python

a

array([[0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [5., 5., 5., 5., 5., 5., 5., 5., 5., 5.],
       [6., 6., 6., 6., 6., 6., 6., 6., 6., 6.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.]])

snippet.python

pkb = np.zeros((3,3))
pkb[1,1] = 1
pkb[0,1] = -1
 
dkb = np.zeros((3,3))
dkb[1,1] = 1
dkb[0,0] = -1/2
 
diffs = np.zeros(shape=(*a.shape, 8))
for rot in range(4):
    pk = np.rot90(pkb, k=rot)
    dk = np.rot90(dkb, k=rot)
    pc = correlate(a, pk, mode='constant', cval=0.0)
    dc = correlate(a, dk, mode='constant', cval=0.0)
    diffs[:, :, rot] = pc
    diffs[:, :, rot+4] = dc

snippet.python

def make_2d_kernels():
    # for the peripheral elements
    pkb = np.zeros((3,3))
    pkb[1,1] = 1
    pkb[0,1] = -1
 
    # for the diagonal elements
    dkb = np.zeros((3,3))
    dkb[1,1] = 1
    dkb[0,0] = -1/2
 
    kernels = np.zeros(shape=(8,3,3))
    for rot in range(4):
        pk = np.rot90(pkb, k=rot)
        dk = np.rot90(dkb, k=rot)
        kernels[rot,:,:] = pk
        kernels[rot+4,:,:] = dk
    return kernels

snippet.python

make_3d_kernels()[2]

array([[[ 0.,  0.,  0.],
        [ 0.,  0.,  0.],
        [ 0.,  0.,  0.]],

       [[ 0., -1.,  0.],
        [ 0.,  1.,  0.],
        [ 0.,  0.,  0.]],

       [[ 0.,  0.,  0.],
        [ 0.,  0.,  0.],
        [ 0.,  0.,  0.]]])

snippet.python

np.rot90(dkb, k=0)

array([[-1.,  0.,  0.],
       [ 0.,  1.,  0.],
       [ 0.,  0.,  0.]])

python