NMM week 7: Discrete Elements¶

Question 1¶

Redo the finite element beam-bending problem (8.2) with a discrete element method, and compare the results. Use a triangular lattice in 2D, with a linear elastic nearest neighbor force law.¶

Instructions: the beam can be a linear spring. Compare this result to the beam-bending problem. In the FEA you get a nice curve. Let it settle to an equilibrium. Do the 2D beam from last week. Additional instructions: play with the force law and add bond damping.

Newton's equations of motion: $$\overrightarrow{F}=m\overrightarrow{a}$$

To solve this, I'm choosing the half-step Verlet. Below are the equations given in the textbook.

Start with the new velocity on the lefthand side, and half Euler step $h\overrightarrow{a}(t)/2$ on the last part of the righthand side: $$\overrightarrow{v}(t+h/2)= \overrightarrow{v}(t)+h\overrightarrow{a}(t)/2$$

New position on the lefthand side: $$\overrightarrow{x}(t+h)= \overrightarrow{x}(t)+h\overrightarrow{v}(t+h/2)$$

Force law on righthand side:
$$\overrightarrow{x}(t+h)= \overrightarrow{a}(t+h)$$

Plug in here, throw out the old $v$, and update the velocity. $$\overrightarrow{v}(t+h)= \overrightarrow{v}(t+h/2)+h\overrightarrow{a}(t+h)/2$$

The way to code this up in Taichi/Python:

  1. Create 3 fields (for position, velocity, force)
  2. A function that initializes positions
  3. A function that calculates forces (Are the particles close? Force attracting vs repelling)
  4. A function that integrates forward in time, updating velocities and positions based on forces.
  5. Visualizing results

Installing Miniconda and Taichi¶

To run Taichi, first install Miniconda. The following steps were taken:

  1. Installing Miniconda.
  2. Make conda environment (name it taichi)
    • To activate conda: conda activate
    • To activate the taichi environment: conda activate taichi
    • To exit: conda deactivate
  3. Installing Taichi here
  4. Git clone Taichi locally

Taichi basics¶

Let's start with the basics of Taichi with this hello world guide: kernels, fields and functions. Here's a powerpoint walkthrough of Taichi.

A kernel is the basic unit for execution in Taichi and the entry point from which Taichi's runtime takes over from Python's virtual machine. You call a kernel the same way as you call a Python function, and you can switch back and forth between Taichi's runtime and Python's virtual machine.

Creating a Taichi field, which is defined as a 'multi-multi-dimensional array of elements, where elements can be: Scalar, Vector, Matrix or Struct. Example of a cloth. GGUI docs with a tutorial on how to create and render 2D and 3D objects.

Project management¶

First spiral: Python + NumPy

  1. Learn NumPy
  2. Visualize in matplotlib or Taichi

NumPy tutorials¶

  • https://www.w3schools.com/python/numpy/default.asp
  • https://www.geeksforgeeks.org/numpy-tutorial/
  • https://realpython.com/numpy-tutorial/
In [ ]:
import taichi as ti
gui = ti.GUI('NMM week 7: beam bending')

# constants
F = 
m = # mass of particle
l = # length of beam
h = # height of beam
n_particles = 128 # grid 

E, nu = 5e3, 0.2  # Young's modulus and Poisson's ratio

# Create 3 arrays (taichi fields - manages data for us in a smart way): 
1 array stores position of particles
1 array stores velocity of all particles
1 array stores forces or accelerations (f=


position_particles = ti.Vector.field(n=2, dtype=float, shape=(n, 2)) # position
velocity_particles = ti.Vector.field(n=2, dtype=float, shape=(n, 2)) # velocity
acceleration_particles = ti.Matrix.field(2, 2, dtype=float, shape=n_particles)  # affine velocity field

# Create initial conditions                           
                                        
# triangular grid = how we initialize positions
Triangular grid: function that generates this grid for us. We want the coordinates of points that would create a triangle. 

# Initialize velocities
We start with all velocities being 0. Taichi is doing this zero by 
Edit data in the 3 taichi fields
                                        
# Visualize results
Taichi rendering code

Python + NumPy¶

Data structures:

  1. Set up triangular grid
  2. Create 3 fields (for position, velocity, force)
  3. A function that initializes positions
  4. A function that calculates forces (Are the particles close? Force attracting vs repelling)
  5. A function that integrates forward in time, updating velocities and positions based on forces.
  6. Visualizing results

Create a triangular grid¶

First, start with setting up the triangular grid of particles that makes up the beam.

Line 20: axes set aspect to equal explaination.

In [36]:
import numpy as np
import matplotlib.pyplot as plt

def triangle_grid_coords(i, j):
    """
    Takes the integers of the points of the first triangle in the grid. 
    """
    tri_x = np.array([1.0, 0.0])
    tri_y = np.array([0.5, np.sqrt(3) /2])
    return i * tri_x + j * tri_y

# Below is a faster way to write a list 
points = np.array([triangle_grid_coords(i,j) for i in range (10) for j in range(5)]) 

print(points.shape) 
print(points)

fig, ax = plt.subplots()
ax.set_aspect("equal") 
ax.scatter(points[:,0], points[:,1], color="pink")
plt.show()

## Access a specific point

i = 1
j = 1
coords = triangle_grid_coords(i, j)
print(f"lattice coordinates:{i},{j}")
print(f"cartesian coordinates:{coords}")

# if you see curly braces, print it as python code 

(50, 2)
[[ 0.          0.        ]
 [ 0.5         0.8660254 ]
 [ 1.          1.73205081]
 [ 1.5         2.59807621]
 [ 2.          3.46410162]
 [ 1.          0.        ]
 [ 1.5         0.8660254 ]
 [ 2.          1.73205081]
 [ 2.5         2.59807621]
 [ 3.          3.46410162]
 [ 2.          0.        ]
 [ 2.5         0.8660254 ]
 [ 3.          1.73205081]
 [ 3.5         2.59807621]
 [ 4.          3.46410162]
 [ 3.          0.        ]
 [ 3.5         0.8660254 ]
 [ 4.          1.73205081]
 [ 4.5         2.59807621]
 [ 5.          3.46410162]
 [ 4.          0.        ]
 [ 4.5         0.8660254 ]
 [ 5.          1.73205081]
 [ 5.5         2.59807621]
 [ 6.          3.46410162]
 [ 5.          0.        ]
 [ 5.5         0.8660254 ]
 [ 6.          1.73205081]
 [ 6.5         2.59807621]
 [ 7.          3.46410162]
 [ 6.          0.        ]
 [ 6.5         0.8660254 ]
 [ 7.          1.73205081]
 [ 7.5         2.59807621]
 [ 8.          3.46410162]
 [ 7.          0.        ]
 [ 7.5         0.8660254 ]
 [ 8.          1.73205081]
 [ 8.5         2.59807621]
 [ 9.          3.46410162]
 [ 8.          0.        ]
 [ 8.5         0.8660254 ]
 [ 9.          1.73205081]
 [ 9.5         2.59807621]
 [10.          3.46410162]
 [ 9.          0.        ]
 [ 9.5         0.8660254 ]
 [10.          1.73205081]
 [10.5         2.59807621]
 [11.          3.46410162]]

lattice coordinates:1,1
cartesian coordinates:[1.5       0.8660254]

That's correct. Now shift the triangular grid to 'fill' the rectangle. Let's first get rid of the code we don't need.

In [35]:
import numpy as np
import matplotlib.pyplot as plt

def triangle_grid_coords(i, j):
    tri_x = np.array([1.0, 0.0])
    tri_y = np.array([0.5, np.sqrt(3) /2])
    return i * tri_x + j * tri_y

points = np.array([triangle_grid_coords(i,j) for i in range (10) for j in range(5)]) 

fig, ax = plt.subplots()
ax.set_aspect("equal") 
ax.scatter(points[:,0], points[:,1], color="pink")
plt.show()

Shifting the coordinates is remarkibly simple, by changing line 9 from:

points = np.array([triangle_grid_coords(i,j) for i in range (10) for j in range(5)])

into:

points = np.array([triangle_grid_coords(i - j // 2, j) for i in range (10) for j in range(5)])

Good to know in Python:

1/2 gives a floating point result.

1 // 2 gives an integer (always rounds down)

In [37]:
import numpy as np
import matplotlib.pyplot as plt

def triangle_grid_coords(i, j):
    tri_x = np.array([1.0, 0.0])
    tri_y = np.array([0.5, np.sqrt(3) /2])
    return i * tri_x + j * tri_y

points = np.array([triangle_grid_coords(i - j // 2, j) for i in range (10) for j in range(5)]) 
# i-j // 2 to shift everything above 1 to the left. 

fig, ax = plt.subplots()
ax.set_aspect("equal") 
ax.scatter(points[:,0], points[:,1], color="pink")
plt.show()

Initial conditions of particles¶

Now we have our beam geometry! Next up: move points in time according to Newton's laws.

  1. What are the forces that act on the points? (internal forces and external forces)
  2. Add Newton's law

In line 28 and 29, I'm using np.zeroes_like for the first time. This is: same shape, only zeros.

In [23]:
import numpy as np
import matplotlib.pyplot as plt

def triangle_grid_coords(i, j):
    tri_x = np.array([1.0, 0.0])
    tri_y = np.array([0.5, np.sqrt(3) /2])
    return i * tri_x + j * tri_y


beam_length = 2    # m
beam_height = 0.3  # m
beam_mass = 10     # kg

n_particles_x = 20
n_particles_y = 10
n_particles = n_particles_x * n_particles_y

particle_mass = beam_mass / n_particles


# Define the initial conditions of the particles 
positions = np.array([
    triangle_grid_coords(i - j // 2, j) 
    for i in range (n_particles_x) 
    for j in range(n_particles_y)
]) 

velocities = np.zeros_like(positions) 
accelerations = np.zeros_like(positions)

# Plot the beam 

fig, ax = plt.subplots()
ax.set_aspect("equal") 
ax.scatter(positions[:,0], positions[:,1], color="pink")
plt.show()

Rewrite lattice coords Indicate where we're going to apply the external force (evenly distributed load)

In [28]:
import numpy as np
import matplotlib.pyplot as plt

def triangle_grid_coords(i, j):
    tri_x = np.array([1.0, 0.0])
    tri_y = np.array([0.5, np.sqrt(3) /2])
    return i * tri_x + j * tri_y

# Define the constants. 
beam_length = 2    # m
beam_height = 0.3  # m
beam_mass = 10     # kg

n_particles_x = 20
n_particles_y = 10
n_particles = n_particles_x * n_particles_y

particle_mass = beam_mass / n_particles



# Define the initial conditions of the simulation - unroll to look at positions and 
lattice_coords = [] # create regular list 
particles_with_external_force = []
particles_without_external_force = []
particle_index = 0 # The row in our position array

for i in range (n_particles_x):
    for j in range(n_particles_y):
        lattice_coords.append( (i - j // 2, j) ) # Add it to the end of the list , and add extra to indicate it's a tuple. 
        
        if i > n_particles_x - 5: # 5 rows of particles the force is applied to 
            particles_with_external_force.append(particle_index)
        else: 
            particles_without_external_force.append(particle_index)
        particle_index += 1
positions = np.array([triangle_grid_coords(i,j) for i, j in lattice_coords])

print(lattice_coords)        
print(positions)


velocities = np.zeros_like(positions) 
accelerations = np.zeros_like(positions)

# Apply an external force to the beam. 




# Plot the beam. 

fig, ax = plt.subplots()
ax.set_aspect("equal") 
ax.scatter(
    positions[particles_without_external_force, 0], 
    positions[particles_without_external_force, 1], 
    color="pink",
)
ax.scatter(
    positions[particles_with_external_force, 0], 
    positions[particles_with_external_force, 1], 
    color="grey",
)
plt.show()

[(0, 0), (0, 1), (-1, 2), (-1, 3), (-2, 4), (-2, 5), (-3, 6), (-3, 7), (-4, 8), (-4, 9), (1, 0), (1, 1), (0, 2), (0, 3), (-1, 4), (-1, 5), (-2, 6), (-2, 7), (-3, 8), (-3, 9), (2, 0), (2, 1), (1, 2), (1, 3), (0, 4), (0, 5), (-1, 6), (-1, 7), (-2, 8), (-2, 9), (3, 0), (3, 1), (2, 2), (2, 3), (1, 4), (1, 5), (0, 6), (0, 7), (-1, 8), (-1, 9), (4, 0), (4, 1), (3, 2), (3, 3), (2, 4), (2, 5), (1, 6), (1, 7), (0, 8), (0, 9), (5, 0), (5, 1), (4, 2), (4, 3), (3, 4), (3, 5), (2, 6), (2, 7), (1, 8), (1, 9), (6, 0), (6, 1), (5, 2), (5, 3), (4, 4), (4, 5), (3, 6), (3, 7), (2, 8), (2, 9), (7, 0), (7, 1), (6, 2), (6, 3), (5, 4), (5, 5), (4, 6), (4, 7), (3, 8), (3, 9), (8, 0), (8, 1), (7, 2), (7, 3), (6, 4), (6, 5), (5, 6), (5, 7), (4, 8), (4, 9), (9, 0), (9, 1), (8, 2), (8, 3), (7, 4), (7, 5), (6, 6), (6, 7), (5, 8), (5, 9), (10, 0), (10, 1), (9, 2), (9, 3), (8, 4), (8, 5), (7, 6), (7, 7), (6, 8), (6, 9), (11, 0), (11, 1), (10, 2), (10, 3), (9, 4), (9, 5), (8, 6), (8, 7), (7, 8), (7, 9), (12, 0), (12, 1), (11, 2), (11, 3), (10, 4), (10, 5), (9, 6), (9, 7), (8, 8), (8, 9), (13, 0), (13, 1), (12, 2), (12, 3), (11, 4), (11, 5), (10, 6), (10, 7), (9, 8), (9, 9), (14, 0), (14, 1), (13, 2), (13, 3), (12, 4), (12, 5), (11, 6), (11, 7), (10, 8), (10, 9), (15, 0), (15, 1), (14, 2), (14, 3), (13, 4), (13, 5), (12, 6), (12, 7), (11, 8), (11, 9), (16, 0), (16, 1), (15, 2), (15, 3), (14, 4), (14, 5), (13, 6), (13, 7), (12, 8), (12, 9), (17, 0), (17, 1), (16, 2), (16, 3), (15, 4), (15, 5), (14, 6), (14, 7), (13, 8), (13, 9), (18, 0), (18, 1), (17, 2), (17, 3), (16, 4), (16, 5), (15, 6), (15, 7), (14, 8), (14, 9), (19, 0), (19, 1), (18, 2), (18, 3), (17, 4), (17, 5), (16, 6), (16, 7), (15, 8), (15, 9)]
[[ 0.          0.        ]
 [ 0.5         0.8660254 ]
 [ 0.          1.73205081]
 [ 0.5         2.59807621]
 [ 0.          3.46410162]
 [ 0.5         4.33012702]
 [ 0.          5.19615242]
 [ 0.5         6.06217783]
 [ 0.          6.92820323]
 [ 0.5         7.79422863]
 [ 1.          0.        ]
 [ 1.5         0.8660254 ]
 [ 1.          1.73205081]
 [ 1.5         2.59807621]
 [ 1.          3.46410162]
 [ 1.5         4.33012702]
 [ 1.          5.19615242]
 [ 1.5         6.06217783]
 [ 1.          6.92820323]
 [ 1.5         7.79422863]
 [ 2.          0.        ]
 [ 2.5         0.8660254 ]
 [ 2.          1.73205081]
 [ 2.5         2.59807621]
 [ 2.          3.46410162]
 [ 2.5         4.33012702]
 [ 2.          5.19615242]
 [ 2.5         6.06217783]
 [ 2.          6.92820323]
 [ 2.5         7.79422863]
 [ 3.          0.        ]
 [ 3.5         0.8660254 ]
 [ 3.          1.73205081]
 [ 3.5         2.59807621]
 [ 3.          3.46410162]
 [ 3.5         4.33012702]
 [ 3.          5.19615242]
 [ 3.5         6.06217783]
 [ 3.          6.92820323]
 [ 3.5         7.79422863]
 [ 4.          0.        ]
 [ 4.5         0.8660254 ]
 [ 4.          1.73205081]
 [ 4.5         2.59807621]
 [ 4.          3.46410162]
 [ 4.5         4.33012702]
 [ 4.          5.19615242]
 [ 4.5         6.06217783]
 [ 4.          6.92820323]
 [ 4.5         7.79422863]
 [ 5.          0.        ]
 [ 5.5         0.8660254 ]
 [ 5.          1.73205081]
 [ 5.5         2.59807621]
 [ 5.          3.46410162]
 [ 5.5         4.33012702]
 [ 5.          5.19615242]
 [ 5.5         6.06217783]
 [ 5.          6.92820323]
 [ 5.5         7.79422863]
 [ 6.          0.        ]
 [ 6.5         0.8660254 ]
 [ 6.          1.73205081]
 [ 6.5         2.59807621]
 [ 6.          3.46410162]
 [ 6.5         4.33012702]
 [ 6.          5.19615242]
 [ 6.5         6.06217783]
 [ 6.          6.92820323]
 [ 6.5         7.79422863]
 [ 7.          0.        ]
 [ 7.5         0.8660254 ]
 [ 7.          1.73205081]
 [ 7.5         2.59807621]
 [ 7.          3.46410162]
 [ 7.5         4.33012702]
 [ 7.          5.19615242]
 [ 7.5         6.06217783]
 [ 7.          6.92820323]
 [ 7.5         7.79422863]
 [ 8.          0.        ]
 [ 8.5         0.8660254 ]
 [ 8.          1.73205081]
 [ 8.5         2.59807621]
 [ 8.          3.46410162]
 [ 8.5         4.33012702]
 [ 8.          5.19615242]
 [ 8.5         6.06217783]
 [ 8.          6.92820323]
 [ 8.5         7.79422863]
 [ 9.          0.        ]
 [ 9.5         0.8660254 ]
 [ 9.          1.73205081]
 [ 9.5         2.59807621]
 [ 9.          3.46410162]
 [ 9.5         4.33012702]
 [ 9.          5.19615242]
 [ 9.5         6.06217783]
 [ 9.          6.92820323]
 [ 9.5         7.79422863]
 [10.          0.        ]
 [10.5         0.8660254 ]
 [10.          1.73205081]
 [10.5         2.59807621]
 [10.          3.46410162]
 [10.5         4.33012702]
 [10.          5.19615242]
 [10.5         6.06217783]
 [10.          6.92820323]
 [10.5         7.79422863]
 [11.          0.        ]
 [11.5         0.8660254 ]
 [11.          1.73205081]
 [11.5         2.59807621]
 [11.          3.46410162]
 [11.5         4.33012702]
 [11.          5.19615242]
 [11.5         6.06217783]
 [11.          6.92820323]
 [11.5         7.79422863]
 [12.          0.        ]
 [12.5         0.8660254 ]
 [12.          1.73205081]
 [12.5         2.59807621]
 [12.          3.46410162]
 [12.5         4.33012702]
 [12.          5.19615242]
 [12.5         6.06217783]
 [12.          6.92820323]
 [12.5         7.79422863]
 [13.          0.        ]
 [13.5         0.8660254 ]
 [13.          1.73205081]
 [13.5         2.59807621]
 [13.          3.46410162]
 [13.5         4.33012702]
 [13.          5.19615242]
 [13.5         6.06217783]
 [13.          6.92820323]
 [13.5         7.79422863]
 [14.          0.        ]
 [14.5         0.8660254 ]
 [14.          1.73205081]
 [14.5         2.59807621]
 [14.          3.46410162]
 [14.5         4.33012702]
 [14.          5.19615242]
 [14.5         6.06217783]
 [14.          6.92820323]
 [14.5         7.79422863]
 [15.          0.        ]
 [15.5         0.8660254 ]
 [15.          1.73205081]
 [15.5         2.59807621]
 [15.          3.46410162]
 [15.5         4.33012702]
 [15.          5.19615242]
 [15.5         6.06217783]
 [15.          6.92820323]
 [15.5         7.79422863]
 [16.          0.        ]
 [16.5         0.8660254 ]
 [16.          1.73205081]
 [16.5         2.59807621]
 [16.          3.46410162]
 [16.5         4.33012702]
 [16.          5.19615242]
 [16.5         6.06217783]
 [16.          6.92820323]
 [16.5         7.79422863]
 [17.          0.        ]
 [17.5         0.8660254 ]
 [17.          1.73205081]
 [17.5         2.59807621]
 [17.          3.46410162]
 [17.5         4.33012702]
 [17.          5.19615242]
 [17.5         6.06217783]
 [17.          6.92820323]
 [17.5         7.79422863]
 [18.          0.        ]
 [18.5         0.8660254 ]
 [18.          1.73205081]
 [18.5         2.59807621]
 [18.          3.46410162]
 [18.5         4.33012702]
 [18.          5.19615242]
 [18.5         6.06217783]
 [18.          6.92820323]
 [18.5         7.79422863]
 [19.          0.        ]
 [19.5         0.8660254 ]
 [19.          1.73205081]
 [19.5         2.59807621]
 [19.          3.46410162]
 [19.5         4.33012702]
 [19.          5.19615242]
 [19.5         6.06217783]
 [19.          6.92820323]
 [19.5         7.79422863]]
In [44]:
import numpy as np
import matplotlib.pyplot as plt

def triangle_grid_coords(i, j):
    tri_x = np.array([1.0, 0.0])
    tri_y = np.array([0.5, np.sqrt(3) /2])
    return i * tri_x + j * tri_y

# Define the beam parameters. 
beam_length = 2    # m
beam_height = 0.3  # m
beam_mass = 10     # kg

# Define the particles. 
lattice_pitch = 0.05 # m
n_particles_x = int(beam_length / lattice_pitch) 
n_particles_y = int(beam_height / (lattice_pitch * np.sqrt(3) /2) ) +1 
n_particles = n_particles_x * n_particles_y
particle_mass = beam_mass / n_particles


# Define the initial conditions of the simulation - unroll to look at positions and 
lattice_coords = [] # create regular list 
particles_with_external_force = []
particles_without_external_force = []
particle_index = 0 # The row in our position array

for i in range (n_particles_x):
    for j in range(n_particles_y):
        lattice_coords.append( (i - j // 2, j) ) # Add it to the end of the list , and add extra to indicate it's a tuple. 
        
        if i > n_particles_x - 5: # 5 rows of particles the force is applied to 
            particles_with_external_force.append(particle_index)
        else: 
            particles_without_external_force.append(particle_index)
        particle_index += 1
positions = lattice_pitch * np.array([triangle_grid_coords(i,j) for i, j in lattice_coords])


velocities = np.zeros_like(positions) 
accelerations = np.zeros_like(positions)


# Plot the beam. 

fig, ax = plt.subplots()
ax.set_aspect("equal") 
ax.scatter(
    positions[particles_without_external_force, 0], 
    positions[particles_without_external_force, 1], 
    color="grey",
)
ax.scatter(
    positions[particles_with_external_force, 0], 
    positions[particles_with_external_force, 1], 
    color="pink",
)
plt.show()

adding forces

In [54]:
import numpy as np
import matplotlib.pyplot as plt

def triangle_grid_coords(i, j):
    tri_x = np.array([1.0, 0.0])
    tri_y = np.array([0.5, np.sqrt(3) /2])
    return i * tri_x + j * tri_y

def apply_external_force(accelerations, particles_with_external_force, force_per_particle, particle_mass):
    for particle_index in particles_with_external_force:
        accelerations[particle_index, 1]  += force_per_particle / particle_mass # 1 to choose column 1 which gives y coordinate   
        
# Define the beam parameters. 
beam_length = 2    # m
beam_height = 0.3  # m
beam_mass = 10     # kg

# Define the particle properties. 
lattice_pitch = 0.05 # m
n_particles_x = int(beam_length / lattice_pitch) 
n_particles_y = int(beam_height / (lattice_pitch * np.sqrt(3) /2) ) +1 
n_particles = n_particles_x * n_particles_y
particle_mass = beam_mass / n_particles

# Define external force. 
total_external_force = -10 # N
force_region_length = 0.1 * beam_length # m 

# Define the initial positions of the particles.
lattice_coords = np.array(
    [(i-j // 2, j) for i in range(n_particles_x) for j in range(n_particles_y)]
)
positions = lattice_pitch * np.array([triangle_grid_coords(i,j) for i, j in lattice_coords])


# Determine which particles are acted on by the external force. 
particles_with_external_force = []
particles_without_external_force = []

for particle_index in range (n_particles):
    j = lattice_coords[particle_index, 1]
    x = positions[particle_index, 0]
    if j == n_particles_y -1 and x > beam_length - force_region_length:
        particles_with_external_force.append(particle_index)
    else: 
        particles_without_external_force.append(particle_index)

# Define initial velocities and accelerations. 
velocities = np.zeros_like(positions) 
accelerations = np.zeros_like(positions)


# Compute the force per particle. 
force_per_particle = total_external_force / len(particles_with_external_force)
print(f"Force per particle: {force_per_particle} N")

apply_external_force(accelerations, particles_with_external_force, force_per_particle, particle_mass)
print(accelerations)

# Plot the beam. 
fig, ax = plt.subplots()
ax.set_aspect("equal") 
ax.scatter(
    positions[particles_without_external_force, 0], 
    positions[particles_without_external_force, 1], 
    color="grey",
)
ax.scatter(
    positions[particles_with_external_force, 0], 
    positions[particles_with_external_force, 1], 
    color="pink",
)
plt.show()

Force per particle: -3.3333333333333335 N
[[  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.        ]
 [  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.        ]
 [  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.        ]
 [  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.        ]
 [  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.        ]
 [  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.        ]
 [  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.        ]
 [  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.        ]
 [  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.        ]
 [  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.        ]
 [  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.         -93.33333333]
 [  0.           0.        ]
 [  0.           0.        ]
 [  0.           0.        ]
 [  0.           0.        ]
 [  0.           0.        ]
 [  0.           0.        ]
 [  0.         -93.33333333]
 [  0.           0.        ]
 [  0.           0.        ]
 [  0.           0.        ]
 [  0.           0.        ]
 [  0.           0.        ]
 [  0.           0.        ]
 [  0.         -93.33333333]]

New function for velocity verlet.

Newton's equations of motion: $$\overrightarrow{F}=m\overrightarrow{a}$$

To solve this, I'm choosing the half-step Verlet equations:

Start with the new velocity on the lefthand side, and half Euler step $h\overrightarrow{a}(t)/2$ on the last part of the righthand side: $$\overrightarrow{v}(t+h/2)= \overrightarrow{v}(t)+h\overrightarrow{a}(t)/2$$

New position on the lefthand side: $$\overrightarrow{x}(t+h)= \overrightarrow{x}(t)+h\overrightarrow{v}(t+h/2)$$

Force law on righthand side:
$$\overrightarrow{x}(t+h)= \overrightarrow{a}(t+h)$$

Plug in here, throw out the old $v$, and update the velocity. $$\overrightarrow{v}(t+h)= \overrightarrow{v}(t+h/2)+h\overrightarrow{a}(t+h)/2$$

In [63]:
import numpy as np
import matplotlib.pyplot as plt

def triangle_grid_coords(i, j):
    tri_x = np.array([1.0, 0.0])
    tri_y = np.array([0.5, np.sqrt(3) /2])
    return i * tri_x + j * tri_y

def apply_external_force(accelerations, particles_with_external_force, force_per_particle, particle_mass): # in parentheses what you're going to call. 
    for particle_index in particles_with_external_force:
        accelerations[particle_index, 1]  += force_per_particle / particle_mass # 1 to choose column 1 which gives y coordinate   

def half_velocity_verlet(
    positions, 
    velocities, 
    accelerations, 
    compute_accelerations, 
    dt,
): 
    # From equation: half-step velocities by h/2 * accelerations.
    velocities = velocities + 0.5 * dt * accelerations
    
    # Update the positions.
    positions = positions + dt * velocities
    
    # Compute new forces. 
    accelerations = compute_accelerations(positions)
    
    # Half-step velocities. 
    velocities = velocities + 0.5 * dt * accelerations
    
    return positions, velocities, accelerations

# Define the beam parameters. 
beam_length = 2    # m
beam_height = 0.3  # m
beam_mass = 10     # kg

# Define the particle properties. 
lattice_pitch = 0.05 # m
n_particles_x = int(beam_length / lattice_pitch) 
n_particles_y = int(beam_height / (lattice_pitch * np.sqrt(3) /2) ) +1 
n_particles = n_particles_x * n_particles_y
particle_mass = beam_mass / n_particles

# Define external force. 
total_external_force = -10 # N
force_region_length = 0.1 * beam_length # m 

# Define time step. 
dt = 0.01 # s 

# Define the initial positions of the particles.
lattice_coords = np.array(
    [(i-j // 2, j) for i in range(n_particles_x) for j in range(n_particles_y)]
)
positions = lattice_pitch * np.array([triangle_grid_coords(i,j) for i, j in lattice_coords])


# Determine which particles are acted on by the external force. 
particles_with_external_force = []
particles_without_external_force = []

for particle_index in range (n_particles):
    j = lattice_coords[particle_index, 1]
    x = positions[particle_index, 0]
    if j == n_particles_y -1 and x > beam_length - force_region_length:
        particles_with_external_force.append(particle_index)
    else: 
        particles_without_external_force.append(particle_index)

# Define initial velocities and accelerations. 
velocities = np.zeros_like(positions) 
accelerations = np.zeros_like(positions)


# Compute the force per particle. 
force_per_particle = total_external_force / len(particles_with_external_force)
print(f"Force per particle: {force_per_particle} N")


def compute_accelerations(positions): 
    accelerations = np.zeros_like(positions)
    apply_external_force(accelerations, particles_with_external_force, force_per_particle, particle_mass)
    print(positions[particles_with_external_force, :])
    return accelerations

# Make 10 steps of velocity verlet. 
for step in range(10):
    positions, velocities, accelerations = half_velocity_verlet(
        positions, 
        velocities, 
        accelerations, 
        compute_accelerations, 
        dt,
    )

# Plot the beam. 
fig, ax = plt.subplots()
ax.set_aspect("equal") 
ax.scatter(
    positions[particles_without_external_force, 0], 
    positions[particles_without_external_force, 1], 
    color="grey",
)
ax.scatter(
    positions[particles_with_external_force, 0], 
    positions[particles_with_external_force, 1], 
    color="pink",
)
plt.show()

Force per particle: -3.3333333333333335 N
[[1.85       0.25980762]
 [1.9        0.25980762]
 [1.95       0.25980762]]
[[1.85       0.25047429]
 [1.9        0.25047429]
 [1.95       0.25047429]]
[[1.85       0.23180762]
 [1.9        0.23180762]
 [1.95       0.23180762]]
[[1.85       0.20380762]
 [1.9        0.20380762]
 [1.95       0.20380762]]
[[1.85       0.16647429]
 [1.9        0.16647429]
 [1.95       0.16647429]]
[[1.85       0.11980762]
 [1.9        0.11980762]
 [1.95       0.11980762]]
[[1.85       0.06380762]
 [1.9        0.06380762]
 [1.95       0.06380762]]
[[ 1.8500000e+00 -1.5257122e-03]
 [ 1.9000000e+00 -1.5257122e-03]
 [ 1.9500000e+00 -1.5257122e-03]]
[[ 1.85       -0.07619238]
 [ 1.9        -0.07619238]
 [ 1.95       -0.07619238]]
[[ 1.85       -0.16019238]
 [ 1.9        -0.16019238]
 [ 1.95       -0.16019238]]

This worked, now let's check if the particles are moving down correctly in each timestep.

In [70]:
import numpy as np
import matplotlib.pyplot as plt

def triangle_grid_coords(i, j):
    tri_x = np.array([1.0, 0.0])
    tri_y = np.array([0.5, np.sqrt(3) /2])
    return i * tri_x + j * tri_y

def apply_external_force(accelerations, particles_with_external_force, force_per_particle, particle_mass): # in parentheses what you're going to call. 
    for particle_index in particles_with_external_force:
        accelerations[particle_index, 1]  += force_per_particle / particle_mass # 1 to choose column 1 which gives y coordinate   

def half_velocity_verlet(
    positions, 
    velocities, 
    accelerations, 
    compute_accelerations, 
    dt,
): 
    # From equation: half-step velocities by h/2 * accelerations.
    velocities = velocities + 0.5 * dt * accelerations
    
    # Update the positions.
    positions = positions + dt * velocities
    
    # Compute new forces. 
    accelerations = compute_accelerations(positions)
    
    # Half-step velocities. 
    velocities = velocities + 0.5 * dt * accelerations
    
    return positions, velocities, accelerations

def plot_beam(positions, particles_with_external_force, particles_without_external_force, step):
    fig, ax = plt.subplots()
    ax.set_aspect("equal") 
    ax.scatter(
        positions[particles_without_external_force, 0], 
        positions[particles_without_external_force, 1], 
        color="grey",
    )
    ax.scatter(
        positions[particles_with_external_force, 0], 
        positions[particles_with_external_force, 1], 
        color="pink",
    )
    plt.show()
    plt.close(fig)


# Define the beam parameters. 
beam_length = 2    # m
beam_height = 0.3  # m
beam_mass = 10     # kg

# Define the particle properties. 
lattice_pitch = 0.05 # m
n_particles_x = int(beam_length / lattice_pitch) 
n_particles_y = int(beam_height / (lattice_pitch * np.sqrt(3) /2) ) +1 
n_particles = n_particles_x * n_particles_y
particle_mass = beam_mass / n_particles

# Define external force. 
total_external_force = -10 # N
force_region_length = 0.1 * beam_length # m 

# Define time step. 
dt = 0.01 # s 

# Define the initial positions of the particles.
lattice_coords = np.array(
    [(i-j // 2, j) for i in range(n_particles_x) for j in range(n_particles_y)]
)
positions = lattice_pitch * np.array([triangle_grid_coords(i,j) for i, j in lattice_coords])


# Determine which particles are acted on by the external force. 
particles_with_external_force = []
particles_without_external_force = []

for particle_index in range (n_particles):
    j = lattice_coords[particle_index, 1]
    x = positions[particle_index, 0]
    if j == n_particles_y -1 and x > beam_length - force_region_length:
        particles_with_external_force.append(particle_index)
    else: 
        particles_without_external_force.append(particle_index)

# Define initial velocities and accelerations. 
velocities = np.zeros_like(positions) 
accelerations = np.zeros_like(positions)


# Compute the force per particle. 
force_per_particle = total_external_force / len(particles_with_external_force)
print(f"Force per particle: {force_per_particle} N")


def compute_accelerations(positions): 
    accelerations = np.zeros_like(positions)
    apply_external_force(accelerations, particles_with_external_force, force_per_particle, particle_mass)
    # print(positions[particles_with_external_force, :])
    return accelerations

# Make 10 steps of velocity verlet. 
for step in range(10):
    positions, velocities, accelerations = half_velocity_verlet(
        positions, 
        velocities, 
        accelerations, 
        compute_accelerations, 
        dt,
    )
    plot_beam(positions, particles_with_external_force, particles_without_external_force, step)
    
    

Force per particle: -3.3333333333333335 N

Adding internal forces¶

Now that the external forces are complete, let's add the internal forces.

$$right force = -k(d-d_0)$$$$left force = k(d-d_0)$$

Where $d_0$ is the lattice_pitch, $d$ is current spring length and $k$ is the spring constant. In 2D, the

$$d = \sqrt{x} $$

Define 6 neighbors around the particle: (i+1, j) (i-1, j) (i, j+1) (i, j-1) (i+1, j-1) (i-1, j+1)

If every particle is responsible for one, so delete 3

In [24]:
import numpy as np
import matplotlib.pyplot as plt

# Define the beam parameters. 
beam_length = 2    # m
beam_height = 0.3  # m
beam_mass = 10     # kg
k = 1000           # N/m

# Define the particle properties. 
lattice_pitch = 0.05 # m
n_particles_x = int(beam_length / lattice_pitch) 
n_particles_y = int(beam_height / (lattice_pitch * np.sqrt(3) /2) ) +1 
n_particles = n_particles_x * n_particles_y
particle_mass = beam_mass / n_particles

# Define external force. 
total_external_force = -1 # N
force_region_length = 0.1 * beam_length # m 

# Define time step. 
dt = 0.002 # s 


def triangle_grid_coords(i, j):
    tri_x = np.array([1.0, 0.0])
    tri_y = np.array([0.5, np.sqrt(3) /2])
    return i * tri_x + j * tri_y

def compute_spring_force(p0, p1, lattice_pitch, k):
    """Compute the force that particle p0 on particle p1."""
    displacement = p1 - p0
    spring_length = np.sqrt(displacement[0] ** 2 + displacement[1] ** 2) # Equation on whiteboard. 
    force_direction = displacement / spring_length
    force_magnitude = -k * (spring_length - lattice_pitch)
    return force_magnitude * force_direction

def apply_internal_force(positions, accelerations, bonds, particle_mass, lattice_pitch, k):
    for particle_index_0, particle_index_1 in bonds: 
        position_0 = positions[particle_index_0]
        position_1 = positions[particle_index_1]
        bond_force = compute_spring_force(position_0, position_1, lattice_pitch, k)
        accelerations[particle_index_0] -= bond_force / particle_mass
        accelerations[particle_index_1] += bond_force / particle_mass

def apply_external_force(accelerations, particles_with_external_force, force_per_particle, particle_mass): # in parentheses what you're going to call. 
    for particle_index in particles_with_external_force:
        accelerations[particle_index, 1]  += force_per_particle / particle_mass # 1 to choose column 1 which gives y coordinate   

def half_velocity_verlet(
    positions, 
    velocities, 
    accelerations, 
    compute_accelerations, 
    dt,
): 
    # From equation: half-step velocities by h/2 * accelerations.
    velocities = velocities + 0.5 * dt * accelerations
    
    # Update the positions.
    positions = positions + dt * velocities
    
    # Compute new forces. 
    accelerations = compute_accelerations(positions)
    
    # Half-step velocities. 
    velocities = velocities + 0.5 * dt * accelerations
    
    return positions, velocities, accelerations

def plot_beam(positions, particles_with_external_force, particles_without_external_force, step):
    fig, ax = plt.subplots()
    ax.set_aspect("equal") 
    ax.scatter(
        positions[particles_without_external_force, 0], 
        positions[particles_without_external_force, 1], 
        color="grey",
    )
    ax.scatter(
        positions[particles_with_external_force, 0], 
        positions[particles_with_external_force, 1], 
        color="pink",
    )
    plt.show()
    plt.close(fig)


# Define the initial positions of the particles.
lattice_coords = np.array(
    [(i-j // 2, j) for i in range(n_particles_x) for j in range(n_particles_y)]
)
positions = lattice_pitch * np.array([triangle_grid_coords(i,j) for i, j in lattice_coords])

# Determine bonds between particles (internal forces). 
lattice_coords_to_particle_index = {(i,j): k for k, (i, j) in enumerate(lattice_coords)} # Creates a 'dictionary'

bonds = []

def add_bond(particle_index, i, j): # Helper function to simplify code. 
    if (i,j) in lattice_coords_to_particle_index:
        other_particle_index = lattice_coords_to_particle_index[(i,j)]
        bonds.append((particle_index, other_particle_index))

for particle_index in range(n_particles):
    i,j = lattice_coords[particle_index]
    add_bond(particle_index, i + 1, j) 
    add_bond(particle_index, i, j + 1)
    add_bond(particle_index, i + 1, j - 1)

# Determine which particles are fixed. 
fixed_particles = []
for particle_index in range(n_particles):
    x = positions[particle_index,0]
    if x <= lattice_pitch:
        fixed_particles.append(particle_index)

# Determine which particles are acted on by the external force. 
particles_with_external_force = []
particles_without_external_force = []

for particle_index in range (n_particles):
    j = lattice_coords[particle_index, 1]
    x = positions[particle_index, 0]
    if j == n_particles_y -1 and x > beam_length - force_region_length:
        particles_with_external_force.append(particle_index)
    else: 
        particles_without_external_force.append(particle_index)

# Define initial velocities and accelerations. 
velocities = np.zeros_like(positions) 
accelerations = np.zeros_like(positions)

# Compute the force per particle. 
force_per_particle = total_external_force / len(particles_with_external_force)
print(f"Force per particle: {force_per_particle} N")


def compute_accelerations(positions): 
    accelerations = np.zeros_like(positions)
    apply_external_force(accelerations, particles_with_external_force, force_per_particle, particle_mass)
    apply_internal_force(positions, accelerations, bonds, particle_mass, lattice_pitch, k)
    accelerations[fixed_particles] = 0.0
    return accelerations

# Make 10 steps of velocity verlet. 
for step in range(100):
    print(f"step{step}")
    for _ in range(10):
        positions, velocities, accelerations = half_velocity_verlet(
            positions, 
            velocities, 
            accelerations, 
            compute_accelerations, 
            dt,
        )
    plot_beam(positions, particles_with_external_force, particles_without_external_force, step)

Force per particle: -0.3333333333333333 N
step0

step1

step2

step3

step4

step5

step6

step7

step8

step9

step10

step11

step12

step13

step14

step15

step16

step17

step18

step19

step20

step21

step22

step23

step24

step25

step26

step27

step28

step29

step30

step31

step32

step33

step34

step35

step36

step37

step38

step39

step40

step41

step42

step43

step44

step45

step46

step47

step48

step49

step50

step51

step52

step53

step54

step55

step56

step57

step58

step59

step60

step61

step62

step63

step64

step65

step66

step67

step68

step69

step70

step71

step72

step73

step74

step75

step76

step77

step78

step79

step80

step81

step82

step83

step84

step85

step86

step87

step88

step89

step90

step91

step92

step93

step94

step95

step96

step97

step98

step99

Question 2¶

Use a discrete element method to simulate filling a container with bodies that are: (a) rigid, (b) elastic, (c) brittle, (d) viscous¶

The goal is to mathematically make a container, then fill it with these 4 materials.

Starting with this Taichi example:

In [ ]:
# Taichi example mpm126.py provided in the library. 

import taichi as ti

ti.init(arch=ti.gpu)  # Try to run on GPU

quality = 1  # Use a larger value for higher-res simulations
n_particles, n_grid = 9000 * quality**2, 128 * quality
dx, inv_dx = 1 / n_grid, float(n_grid)
dt = 1e-4 / quality
p_vol, p_rho = (dx * 0.5)**2, 1
p_mass = p_vol * p_rho
E, nu = 5e3, 0.2  # Young's modulus and Poisson's ratio
mu_0, lambda_0 = E / (2 * (1 + nu)), E * nu / (
    (1 + nu) * (1 - 2 * nu))  # Lame parameters

x = ti.Vector.field(2, dtype=float, shape=n_particles)  # position
v = ti.Vector.field(2, dtype=float, shape=n_particles)  # velocity
C = ti.Matrix.field(2, 2, dtype=float,
                    shape=n_particles)  # affine velocity field
F = ti.Matrix.field(2, 2, dtype=float,
                    shape=n_particles)  # deformation gradient
material = ti.field(dtype=int, shape=n_particles)  # material id
Jp = ti.field(dtype=float, shape=n_particles)  # plastic deformation
grid_v = ti.Vector.field(2, dtype=float,
                         shape=(n_grid, n_grid))  # grid node momentum/velocity
grid_m = ti.field(dtype=float, shape=(n_grid, n_grid))  # grid node mass
gravity = ti.Vector.field(2, dtype=float, shape=())
attractor_strength = ti.field(dtype=float, shape=())
attractor_pos = ti.Vector.field(2, dtype=float, shape=())


@ti.kernel
def substep():
    for i, j in grid_m:
        grid_v[i, j] = [0, 0]
        grid_m[i, j] = 0
    for p in x:  # Particle state update and scatter to grid (P2G)
        base = (x[p] * inv_dx - 0.5).cast(int)
        fx = x[p] * inv_dx - base.cast(float)
        # Quadratic kernels  [http://mpm.graphics   Eqn. 123, with x=fx, fx-1,fx-2]
        w = [0.5 * (1.5 - fx)**2, 0.75 - (fx - 1)**2, 0.5 * (fx - 0.5)**2]
        # deformation gradient update
        F[p] = (ti.Matrix.identity(float, 2) + dt * C[p]) @ F[p]
        # Hardening coefficient: snow gets harder when compressed
        h = ti.max(0.1, ti.min(5, ti.exp(10 * (1.0 - Jp[p]))))
        if material[p] == 1:  # jelly, make it softer
            h = 0.3
        mu, la = mu_0 * h, lambda_0 * h
        if material[p] == 0:  # liquid
            mu = 0.0
        U, sig, V = ti.svd(F[p])
        J = 1.0
        for d in ti.static(range(2)):
            new_sig = sig[d, d]
            if material[p] == 2:  # Snow
                new_sig = min(max(sig[d, d], 1 - 2.5e-2),
                              1 + 4.5e-3)  # Plasticity
            Jp[p] *= sig[d, d] / new_sig
            sig[d, d] = new_sig
            J *= new_sig
        if material[p] == 0:
            # Reset deformation gradient to avoid numerical instability
            F[p] = ti.Matrix.identity(float, 2) * ti.sqrt(J)
        elif material[p] == 2:
            # Reconstruct elastic deformation gradient after plasticity
            F[p] = U @ sig @ V.transpose()
        stress = 2 * mu * (F[p] - U @ V.transpose()) @ F[p].transpose(
        ) + ti.Matrix.identity(float, 2) * la * J * (J - 1)
        stress = (-dt * p_vol * 4 * inv_dx * inv_dx) * stress
        affine = stress + p_mass * C[p]
        for i, j in ti.static(ti.ndrange(3, 3)):
            # Loop over 3x3 grid node neighborhood
            offset = ti.Vector([i, j])
            dpos = (offset.cast(float) - fx) * dx
            weight = w[i][0] * w[j][1]
            grid_v[base + offset] += weight * (p_mass * v[p] + affine @ dpos)
            grid_m[base + offset] += weight * p_mass
    for i, j in grid_m:
        if grid_m[i, j] > 0:  # No need for epsilon here
            # Momentum to velocity
            grid_v[i, j] = (1 / grid_m[i, j]) * grid_v[i, j]
            grid_v[i, j] += dt * gravity[None] * 30  # gravity
            dist = attractor_pos[None] - dx * ti.Vector([i, j])
            grid_v[i, j] += \
                dist / (0.01 + dist.norm()) * attractor_strength[None] * dt * 100
            if i < 3 and grid_v[i, j][0] < 0:
                grid_v[i, j][0] = 0  # Boundary conditions
            if i > n_grid - 3 and grid_v[i, j][0] > 0:
                grid_v[i, j][0] = 0
            if j < 3 and grid_v[i, j][1] < 0:
                grid_v[i, j][1] = 0
            if j > n_grid - 3 and grid_v[i, j][1] > 0:
                grid_v[i, j][1] = 0
    for p in x:  # grid to particle (G2P)
        base = (x[p] * inv_dx - 0.5).cast(int)
        fx = x[p] * inv_dx - base.cast(float)
        w = [0.5 * (1.5 - fx)**2, 0.75 - (fx - 1.0)**2, 0.5 * (fx - 0.5)**2]
        new_v = ti.Vector.zero(float, 2)
        new_C = ti.Matrix.zero(float, 2, 2)
        for i, j in ti.static(ti.ndrange(3, 3)):
            # loop over 3x3 grid node neighborhood
            dpos = ti.Vector([i, j]).cast(float) - fx
            g_v = grid_v[base + ti.Vector([i, j])]
            weight = w[i][0] * w[j][1]
            new_v += weight * g_v
            new_C += 4 * inv_dx * weight * g_v.outer_product(dpos)
        v[p], C[p] = new_v, new_C
        x[p] += dt * v[p]  # advection


@ti.kernel
def reset():
    group_size = n_particles // 3
    for i in range(n_particles):
        x[i] = [
            ti.random() * 0.2 + 0.3 + 0.10 * (i // group_size),
            ti.random() * 0.2 + 0.05 + 0.32 * (i // group_size)
        ]
        material[i] = i // group_size  # 0: fluid 1: jelly 2: snow
        v[i] = [0, 0]
        F[i] = ti.Matrix([[1, 0], [0, 1]])
        Jp[i] = 1
        C[i] = ti.Matrix.zero(float, 2, 2)


print(
    "[Hint] Use WSAD/arrow keys to control gravity. Use left/right mouse buttons to attract/repel. Press R to reset."
)
gui = ti.GUI("Taichi MLS-MPM-128", res=512, background_color=0x112F41)
reset()
gravity[None] = [0, -1]

for frame in range(20000):
    if gui.get_event(ti.GUI.PRESS):
        if gui.event.key == 'r':
            reset()
        elif gui.event.key in [ti.GUI.ESCAPE, ti.GUI.EXIT]:
            break
    if gui.event is not None:
        gravity[None] = [0, 0]  # if had any event
    if gui.is_pressed(ti.GUI.LEFT, 'a'):
        gravity[None][0] = -1
    if gui.is_pressed(ti.GUI.RIGHT, 'd'):
        gravity[None][0] = 1
    if gui.is_pressed(ti.GUI.UP, 'w'):
        gravity[None][1] = 1
    if gui.is_pressed(ti.GUI.DOWN, 's'):
        gravity[None][1] = -1
    mouse = gui.get_cursor_pos()
    gui.circle((mouse[0], mouse[1]), color=0x336699, radius=15)
    attractor_pos[None] = [mouse[0], mouse[1]]
    attractor_strength[None] = 0
    if gui.is_pressed(ti.GUI.LMB):
        attractor_strength[None] = 1
    if gui.is_pressed(ti.GUI.RMB):
        attractor_strength[None] = -1
    for s in range(int(2e-3 // dt)):
        substep()
    gui.circles(x.to_numpy(),
                radius=1.5,
                palette=[0x068587, 0xED553B, 0xEEEEF0],
                palette_indices=material)

    # Change to gui.show(f'{frame:06d}.png') to write images to disk
    gui.show()

For multiple boxes¶

Compared to the first question, the code can be almost the same, but there are no more bonds, all forces need to go to 0, and the shape of the force law depends on the material (e.g. liquid is nonlinear). Pay attention to:

  • Apply internal force (similar to first assignment)
  • Set a cutoff distance for each particle. If distance > cutoff distance the particles will no longer be considered for the interaction.
  • Stability limit depends on the material.

Click here to go back to the main page.