Final Project: Discrete elements for simulating failure in voxel structures

The Nature of Mathematical Modeling

May 22, 2023

SEE https://docs.google.com/presentation/d/10FvxsZxgiPSxUl4yleZDDsTmKQ76YofL3bJn4FKs76I/edit?usp=sharing FOR PRESENTATION

alternately, see: https://fab.cba.mit.edu/classes/864.23/people/Miana/NMM_final.pdf for the PDF version, though note that the presentation contains animations which are not viewable in the pdf.

The code below is scratch work.

All other code and relevant files can be found here: https://gitlab.cba.mit.edu/classes/864.23/site/-/tree/main/people/Miana/py/final This includes python code as well as stl test files.

The summary slide is copied below for reference, but please refer to the presentation link for documentation:

What does it do? Discrete element simulation of voxels, aimed at examining failure modes

Who's done what beforehand? There are a lot of tools for simulating these structures, but not so many focused specifically on failure. The DEM work from Erik and Neil in the DEM chapter is the most direct relationship.

What materials and components were used? This is written in Python and makes use of taichi for acceleration, pyvista for handling meshes, numpy/scipy for various array calculations, and matplotlib for visualization.

What processes were used? This mostly draws on the discrete elements chapter, though that chapter in turn does draw on the ones before it.

What questions were answered? A lot of my personal questions about setting up a discrete element method were answered. Working toward an answer for “is DEM suitable for simulated buckling/fracture/etc in voxel structures?”

How was it evaluated? Mostly qualitatively currently, comparing against real measured results

What are the implications? If it worked this could be a useful prototyping tool for voxel design.

Scratch work v

import pyvista as pv
import numpy as np
from scipy.spatial import KDTree
import matplotlib.pyplot as plt

mesh = pv.read('final/dogbone_type_5.stl')

voxels = pv.voxelize(mesh, density=2.0) # converts to voxels
# 2 gives us 116 points

glyphs = voxels.glyph(factor=0.75, geom=pv.Sphere()) # this converts the geometry to other primitives, where factor scales it.

centers = voxels.cell_centers().points

N = centers.shape[0] # number of particles

# system parameters
total_mass = 0.5 # kg
gravity = 10 # m/s
diss = 1e-4 # dissipation

# particle parameters
p = np.copy(centers) # particle positions <- this is what we'll update
v = np.zeros_like(p) # velocities
a = np.zeros_like(p) # accelerations

# material parameters
mass = total_mass/N # even masses
m = mass*np.ones_like(p)
k = 1e8 # N/m spring constant... maybe modulus actually

# simulation parameters
dt = 0.0001 # time step

print(centers[0])

[1. 1. 1.]

dist_tree = KDTree(p)

dd2,ii2 = dist_tree.query(p,k=2) #get the nearest neighbor, dd2 will be [distance to self (0), and lattice pitch]
pitch = dd2[0,1]
eps = 0.001
q_radius = np.sqrt(2)*pitch + eps # radius we should be within

ind = dist_tree.query_ball_tree(dist_tree, q_radius)
# print(ind)
# print(len(ind))


dd, ii = dist_tree.query(p,k=9,distance_upper_bound=q_radius)
# print(dd)
# print(ii)


# print(dd)
# print(dd2.shape)
# print(dd2)
# print(ii)


# fig, ax = plt.subplots()
# ax.plot(centers[:,0],centers[:,1],'ro')
# ax.plot(centers[ii[0],0],centers[ii[0],1],'ko',alpha = 0.5)
# ax.axis('equal')
# plt.show()

def get3Distance(p,ind):
    N = len(p)
    D = []
    for n in range(N):
        i = ind[n] # list of nearest neighbor's indeces
        amt = len(i)
        dists = np.zeros((amt,3)) # x,y,z for each component's distances
        particles = p[i]
        #print(particles)
        for j in range(amt):
            dists[j] = particles[j] - particles[0]
        D.append(dists)
    return D
init_dist = get3Distance(p,ind)

def force_1D(d, di,bounds):
    force = 0
    distance = d - di
    if di==0:
        force = 0
    else:
        if abs(distance/di) <= bounds[0]:
            force = -k*distance
        elif abs(distance/di) <= bounds[1]:
            #spring_max = -k*bounds[0]*di
            if distance < 0:
                force = k*bounds[0]*di
            else:
                force = -k*bounds[0]*di
#         elif abs(distance/di) <= bounds[2]:
#             pass
        elif distance/di <= -bounds[1]:
            force = k*bounds[0]*di
    #     if abs(distance/di) <= 0.85:
    #         force = 1e7 #something
    #     elif abs(distance/di) < 1.15:
    #         #apply spring force
    #         spring = -k*distance
    #         force = spring
    return force
bounds = [0.15,0.25,0.35]
dist = np.linspace(0,2,100)
forces = np.zeros_like(dist)
for i in range(len(forces)):
    forces[i] = force_1D(dist[i],1,bounds)
# forces = force_1D(dist,di)
plt.plot((dist-1),forces)
plt.xlabel("fractional distance from neutral")
plt.ylabel("force")
plt.show()

def forceLaw(d, di, force, grav=False, dis=1e-5):
    """ 
    force law per particle, as a function of:
    d - distance current
    di - initial/neutral distance
    force - force law
    dis - dissipation
    returns:
    f - force on particle from nearest neighbors (array of 3-vectors of length neighbor amount)
    """
    f = np.zeros_like(d)
    N = d.shape[0]
    dim = d.shape[1]
    bounds = [0.15,0.25,0.35]
    for i in range(N):
        for j in range(dim):
            f[i][j] = force(d[i][j],di[i][j],bounds)
            if grav:
                f[i][j] += gravity
    f = np.sum(f,axis=0)
    return f
force_law = forceLaw(init_dist[0],init_dist[0],force_1D,grav=False)
print(force_law)

[0. 0. 0.]

def integrateForward(p,v,a,mass,init_dist,forceLaw,dt,dis=1e-5):
    di = init_dist
    d = get3Distance(p,ind)
    N = len(d)
    for n in range(N):
        #print(1/m*forceLaw(d[n],di[n],force_1D,grav=False))
        a[n] = 1/mass*forceLaw(d[n],di[n],force_1D,grav=False)
    a += -dis*a #dissipation
    v += a*dt
    p += v*dt
    return p,v,a
p,v,a = integrateForward(p,v,a,mass,init_dist,forceLaw,dt)

tai chi scratch

import pyvista as pv
import numpy as np
from scipy.spatial import KDTree
import matplotlib.pyplot as plt
import matplotlib.animation as animation
import taichi as ti

[Taichi] version 1.4.1, llvm 15.0.5, commit e67c674e, osx, python 3.10.10
[I 05/19/23 15:20:27.723 4577298] [shell.py:_shell_pop_print@23] Graphical python shell detected, using wrapped sys.stdout

ti.init(arch=ti.cpu) # initialize taichi


# IMPORT GEOMETRY
mesh = pv.read('final/dogbone_type_5.stl')
print("File import complete")

# CONVERT TO PARTICLES
voxels = pv.voxelize(mesh, density=1.5) # converts to voxels
# 2 gives us 116 points
print("Particle conversion complete")

glyphs = voxels.glyph(factor=0.75, geom=pv.Sphere()) # this converts the geometry to other primitives, where factor scales it.

centers = voxels.cell_centers().points

N = centers.shape[0] # number of particles

print(f"Particle amount: {N}")

[Taichi] Starting on arch=x64
File import complete
Particle conversion complete
Particle amount: 448

# brute force find the nearest neighbors to instantiate system.

dist_tree = KDTree(centers)
dd2,ii2 = dist_tree.query(centers,k=2) #get the nearest neighbor, dd2 will be [distance to self (0), and lattice pitch]
pitch = dd2[0,1]
eps = 0.001
q_radius = np.sqrt(2)*pitch + eps # radius we should be within
# print(q_radius)

dd, ii = dist_tree.query(centers,k=27,distance_upper_bound=q_radius)
# print(dd)
# print(ii[np.nonzero(dd==np.inf)])
# print(ii)
print(ii[0])
print(centers[ii[0][:7]])
print(centers.shape)

[  0   1   6 224   7 225 230 448 448 448 448 448 448 448 448 448 448 448
 448 448 448 448 448 448 448 448 448]
[[0.75 0.75 0.75]
 [0.75 2.25 0.75]
 [2.25 0.75 0.75]
 [0.75 0.75 2.25]
 [2.25 2.25 0.75]
 [0.75 2.25 2.25]
 [2.25 0.75 2.25]]
(448, 3)

# system parameters
total_mass = 0.5 # kg
gravity = 10 # m/s
diss = 1e-4 # dissipation

# material parameters
mass = total_mass/N # even masses
k = 1e2 # N/m spring constant... maybe modulus actually

# simulation parameters
dt = 0.001 # time step


# taichi vector fields for position, velocity, acceleration
p = ti.Vector.field(n=3,dtype=float, shape=(N,)) # positions
p.from_numpy(centers) # population position with values from centers

v = ti.Vector.field(n=3,dtype=float, shape=(N,)) # velocities
a = ti.Vector.field(n=3,dtype=float, shape=(N,)) # accelerations
f = ti.Vector.field(n=3,dtype=float, shape=(N,)) # forces

# distance tree
dist_tree = KDTree(centers)

dd,ii = dist_tree.query(centers,k=9) # nearest neighbors in 3D

ind = ti.Vector.field(n=9,dtype=int,shape=(N,))
ind.from_numpy(ii)

d = ti.Vector.field(n=3,dtype=float,shape=(N,9))

di = ti.Vector.field(n=3,dtype=float,shape=(N,9)) # this is the initial or neutral distance between points

x_min = np.amin(centers,axis=0)[0]
# print(x_min)
x_max = np.amax(centers,axis=0)[0]
# print(x_max)
edge_index = np.nonzero(centers[:,0]==x_max)[0]
# print(edge_index)
ei = ti.field(ti.i32,edge_index.shape)
ei.from_numpy(edge_index)

@ti.kernel
def getNeutralDistance():
    dim = di.shape[1]
    for n in range(N):
        for i in range(dim):
            di[n,i] = p[ind[n][i]] - p[ind[n][0]]

getNeutralDistance() # initialize the neutral distances

@ti.func
def getDistance():
    dim = d.shape[1]
    for n in range(N):
        for i in range(dim):
            d[n,i] = (p[ind[n][i]] - p[ind[n][0]])-di[n,i]
            
@ti.func
def getSpringForce():
    dim = d.shape[1]
    for i in f:
        for j in range(dim):
            f[i] -= -k*d[i,j] + 0.5*v[i]

@ti.kernel
def stepForward():
    # get forces
    # spring force
    # update RELATIVE displacements d
    getDistance()
    getSpringForce()
    for i in ei:
        f[ei[i]] += [10,0,0]
    for i in f:
        # gravity
        # f[i] -= [0,gravity,0] # gravity in y
        # dissipation
        f[i] -= diss*f[i]
    # apply loads
    # integrate forward
    for i in p:
        a[i] = 1/mass * f[i]
        v[i] += a[i]*dt
        # boundary conditions
        if p[i][0] == x_min:
            p[i] = p[i]
        else:
            p[i] += v[i]*dt

# system parameters
total_mass = 0.5 # kg
gravity = 10 # m/s
diss = 1e-4 # dissipation

# material parameters
mass = total_mass/N # even masses
k = 1e2 # N/m spring constant... maybe modulus actually

# simulation parameters
dt = 0.001 # time step


# taichi vector fields for position, velocity, acceleration
p = ti.Vector.field(n=3,dtype=float, shape=(N,)) # positions
p.from_numpy(centers) # population position with values from centers

v = ti.Vector.field(n=3,dtype=float, shape=(N,)) # velocities
a = ti.Vector.field(n=3,dtype=float, shape=(N,)) # accelerations
f = ti.Vector.field(n=3,dtype=float, shape=(N,)) # forces

# distance tree
dist_tree = KDTree(centers)

dd,ii = dist_tree.query(centers,k=9) # nearest neighbors in 3D

ind = ti.Vector.field(n=9,dtype=int,shape=(N,))
ind.from_numpy(ii)

d = ti.Vector.field(n=3,dtype=float,shape=(N,9))

di = ti.Vector.field(n=3,dtype=float,shape=(N,9)) # this is the initial or neutral distance between points

x_min = np.amin(centers,axis=0)[0]
# print(x_min)
x_max = np.amax(centers,axis=0)[0]
# print(x_max)
edge_index = np.nonzero(centers[:,0]==x_max)[0]
# print(edge_index)
ei = ti.field(ti.i32,edge_index.shape)
ei.from_numpy(edge_index)

@ti.kernel
def getNeutralDistance():
    dim = di.shape[1]
    for n in range(N):
        for i in range(dim):
            di[n,i] = p[ind[n][i]] - p[ind[n][0]]

getNeutralDistance() # initialize the neutral distances

@ti.func
def getDistance():
    dim = d.shape[1]
    for n in range(N):
        for i in range(dim):
            d[n,i] = (p[ind[n][i]] - p[ind[n][0]])-di[n,i]
            
@ti.func
def getSpringForce():
    dim = d.shape[1]
    for i in f:
        for j in range(dim):
            f[i] -= -k*d[i,j]# + 0.5*v[i]

@ti.kernel
def stepForward():
    # get forces
    # spring force
    # update RELATIVE displacements d
    getDistance()
    getSpringForce()
    for i in ei:
        f[ei[i]] += [0.01,0,0]
    for i in f:
        # gravity
        # f[i] -= [0,gravity,0] # gravity in y
        # dissipation
        pass
        # f[i] -= diss*f[i]
    # apply loads
    # integrate forward
    for i in p:
        a[i] = 1/mass * f[i]
        v[i] += a[i]*dt
        # boundary conditions
        if p[i][0] == x_min:
            p[i] = p[i]
        else:
            p[i] += v[i]*dt
            
fig, ax = plt.subplots(subplot_kw=dict(projection='3d'))
# print(p)
plot_p = p.to_numpy()
# print(plot_p)
ax.plot(plot_p[:,0],plot_p[:,1],plot_p[:,2],'k.')
ax.axis('equal')
for i in range(30):
    stepForward()
# print(p)
plot_p = p.to_numpy()
ax.plot(plot_p[:,0],plot_p[:,1],plot_p[:,2],'r.')

plt.show()

[W 05/19/23 15:30:20.048 4577298] [type_check.cpp:type_check_store@36] [$31660] Global store may lose precision: f32 <- f64
File "/Users/mianasmith/Documents/2_NMM/ballz/venv/lib/python3.10/site-packages/taichi/_kernels.py", line 189, in ext_arr_to_matrix:
                        mat[I][p] = arr[I, p]
                        ^^^^^^^^^^^^^^^^^^^^^

[W 05/19/23 15:30:20.048 4577298] [type_check.cpp:type_check_store@36] [$31683] Global store may lose precision: f32 <- f64
File "/Users/mianasmith/Documents/2_NMM/ballz/venv/lib/python3.10/site-packages/taichi/_kernels.py", line 189, in ext_arr_to_matrix:
                        mat[I][p] = arr[I, p]
                        ^^^^^^^^^^^^^^^^^^^^^

[W 05/19/23 15:30:20.048 4577298] [type_check.cpp:type_check_store@36] [$31706] Global store may lose precision: f32 <- f64
File "/Users/mianasmith/Documents/2_NMM/ballz/venv/lib/python3.10/site-packages/taichi/_kernels.py", line 189, in ext_arr_to_matrix:
                        mat[I][p] = arr[I, p]
                        ^^^^^^^^^^^^^^^^^^^^^

[W 05/19/23 15:30:20.139 4577298] [type_check.cpp:type_check_store@36] [$31925] Global store may lose precision: i32 <- i64
File "/Users/mianasmith/Documents/2_NMM/ballz/venv/lib/python3.10/site-packages/taichi/_kernels.py", line 189, in ext_arr_to_matrix:
                        mat[I][p] = arr[I, p]
                        ^^^^^^^^^^^^^^^^^^^^^

[W 05/19/23 15:30:20.139 4577298] [type_check.cpp:type_check_store@36] [$31948] Global store may lose precision: i32 <- i64
File "/Users/mianasmith/Documents/2_NMM/ballz/venv/lib/python3.10/site-packages/taichi/_kernels.py", line 189, in ext_arr_to_matrix:
                        mat[I][p] = arr[I, p]
                        ^^^^^^^^^^^^^^^^^^^^^

[W 05/19/23 15:30:20.139 4577298] [type_check.cpp:type_check_store@36] [$31971] Global store may lose precision: i32 <- i64
File "/Users/mianasmith/Documents/2_NMM/ballz/venv/lib/python3.10/site-packages/taichi/_kernels.py", line 189, in ext_arr_to_matrix:
                        mat[I][p] = arr[I, p]
                        ^^^^^^^^^^^^^^^^^^^^^

[W 05/19/23 15:30:20.139 4577298] [type_check.cpp:type_check_store@36] [$31994] Global store may lose precision: i32 <- i64
File "/Users/mianasmith/Documents/2_NMM/ballz/venv/lib/python3.10/site-packages/taichi/_kernels.py", line 189, in ext_arr_to_matrix:
                        mat[I][p] = arr[I, p]
                        ^^^^^^^^^^^^^^^^^^^^^

[W 05/19/23 15:30:20.139 4577298] [type_check.cpp:type_check_store@36] [$32017] Global store may lose precision: i32 <- i64
File "/Users/mianasmith/Documents/2_NMM/ballz/venv/lib/python3.10/site-packages/taichi/_kernels.py", line 189, in ext_arr_to_matrix:
                        mat[I][p] = arr[I, p]
                        ^^^^^^^^^^^^^^^^^^^^^

[W 05/19/23 15:30:20.139 4577298] [type_check.cpp:type_check_store@36] [$32040] Global store may lose precision: i32 <- i64
File "/Users/mianasmith/Documents/2_NMM/ballz/venv/lib/python3.10/site-packages/taichi/_kernels.py", line 189, in ext_arr_to_matrix:
                        mat[I][p] = arr[I, p]
                        ^^^^^^^^^^^^^^^^^^^^^

[W 05/19/23 15:30:20.139 4577298] [type_check.cpp:type_check_store@36] [$32063] Global store may lose precision: i32 <- i64
File "/Users/mianasmith/Documents/2_NMM/ballz/venv/lib/python3.10/site-packages/taichi/_kernels.py", line 189, in ext_arr_to_matrix:
                        mat[I][p] = arr[I, p]
                        ^^^^^^^^^^^^^^^^^^^^^

[W 05/19/23 15:30:20.139 4577298] [type_check.cpp:type_check_store@36] [$32086] Global store may lose precision: i32 <- i64
File "/Users/mianasmith/Documents/2_NMM/ballz/venv/lib/python3.10/site-packages/taichi/_kernels.py", line 189, in ext_arr_to_matrix:
                        mat[I][p] = arr[I, p]
                        ^^^^^^^^^^^^^^^^^^^^^

[W 05/19/23 15:30:20.139 4577298] [type_check.cpp:type_check_store@36] [$32109] Global store may lose precision: i32 <- i64
File "/Users/mianasmith/Documents/2_NMM/ballz/venv/lib/python3.10/site-packages/taichi/_kernels.py", line 189, in ext_arr_to_matrix:
                        mat[I][p] = arr[I, p]
                        ^^^^^^^^^^^^^^^^^^^^^

[W 05/19/23 15:30:20.214 4577298] [type_check.cpp:type_check_store@36] [$32491] Global store may lose precision: i32 <- i64
File "/Users/mianasmith/Documents/2_NMM/ballz/venv/lib/python3.10/site-packages/taichi/_kernels.py", line 126, in ext_arr_to_tensor:
        tensor[I] = arr[I]
        ^^^^^^^^^^^^^^^^^^

uh.... no

# don't actually think we want this
def getNeutralDistance(p,ind,dd):
    N = len(p)
    dists = []
    for i in range(N):
        index = len(ind[i])
        dists.append(dd[i,:index])
    return dists
NeutralDistance = getNeutralDistance(p, ind,dd)

def getForces(p,NeutralDistance,ind, force_law):
    N = len(p)
    F = np.zeros((N,3)) # forces in direction
    # apply 
force_law = 0
getForces(p,NeutralDistance,ind,force_law)