Springs Mechanics
In this notebook, we'll be building a simple simulation of points connected through springs (inspired by this paper).
Let's first define a triangle of points, as well as springs that connect the points.
In [1]:
import numpy as np
x = np.array([[0., 0.],
[1., 0.],
[0, -1.]])
spring_anchors = np.array([[0, 1],
[1, 2],
[2, 0]])
In [2]:
import holoviews as hv
hv.extension('bokeh', 'matplotlib')
def plot_points(x, spring_anchors):
"""Plots all x as 2D points, connected by springs indicated by spring_anchors."""
lines = []
for i in range(spring_anchors.shape[0]):
a, b = spring_anchors[i]
lines.append(hv.Curve(x[[a, b], :]))
lines = hv.Overlay(lines)
scatter = hv.Scatter((x[:, 0], x[:, 1])).opts(hv.opts.Scatter(size=30, padding=0.2))
labels = hv.Labels({('x', 'y'): x, 'text': [str(n) for n in range(x.shape[0])]}, ['x', 'y'], 'text')
return scatter * lines * labels
plot_points(x, spring_anchors)
Out[2]:
Now, we're going to define springs parameters and masses.
In [3]:
spring_stiffnesses = [2., 2., 2.]
spring_lengths = [.5, .5, .5]
masses = np.ones_like(spring_stiffnesses)
To move the particles, we need to add the spring forces to a force vector and step it forward in time.
Let's first do the force computation.
In [4]:
def compute_spring_forces(x, spring_anchors, spring_lengths, spring_stiffnesses):
"""Computes the spring forces."""
forces = np.zeros_like(x)
for i in range(spring_anchors.shape[0]):
a, b = spring_anchors[i]
x_a, x_b = x[a], x[b]
dist = x_a - x_b
length = np.linalg.norm(dist) + 1e-4
F = (length - spring_lengths[i]) * spring_stiffnesses[i] * dist / length
forces[a] += -F
forces[b] += F
return forces
forces = compute_spring_forces(x, spring_anchors, spring_lengths, spring_stiffnesses)
forces
Out[4]:
array([[ 1.00009999, -1.00009999],
[-2.29304321, -1.29294322],
[ 1.29294322, 2.29304321]])
Now, let's write the time integrator:
In [5]:
def time_step_integration(x, v, forces, masses, dt, damping):
"""Integrates the forces using semi-implicit Euler time integration with damping."""
s = np.exp(-dt * damping)
v_new = s * v + dt * forces / masses.reshape(-1, 1)
x_new = x + dt * v_new
return x_new, v_new
v = np.zeros_like(x)
time_step_integration(x, v, forces, masses, dt=0.1, damping=0.1)
Out[5]:
(array([[ 0.010001 , -0.010001 ],
[ 0.97706957, -0.01292943],
[ 0.01292943, -0.97706957]]), array([[ 0.10001 , -0.10001 ],
[-0.22930432, -0.12929432],
[ 0.12929432, 0.22930432]]))
We can write the forward program:
In [6]:
def forward(x, steps=200, dt=0.5, damping=0.1):
new_x = x.copy()
new_v = np.zeros_like(x)
snapshots = []
snapshots.append(new_x.copy())
for t in range(steps):
forces = compute_spring_forces(new_x, spring_anchors, spring_lengths, spring_stiffnesses)
new_x, new_v = time_step_integration(new_x, new_v, forces, masses, dt=dt, damping=damping)
snapshots.append(new_x.copy())
return snapshots
snapshots = forward(x)
plot_points(snapshots[-1], spring_anchors)
Out[6]:
We can also show some of the coordinate plots over time.
In [7]:
def make_time_plots(snapshots, node_number):
return hv.Curve([snapshot[node_number, 0] for snapshot in snapshots],
label='node{}_x'.format(node_number), kdims='time', vdims='coordinate') * \
hv.Curve([snapshot[node_number, 1] for snapshot in snapshots],
label='node{}_y'.format(node_number), kdims='time', vdims='coordinate')
make_time_plots(snapshots, 0)
Out[7]:
With all that done, we can now do an animation showing the spring evolution over time.
In [8]:
def make_animation(snapshots, n_frames):
hmap = hv.HoloMap({i: plot_points(snapshots[i], spring_anchors) for i in range(0, len(snapshots), n_frames)},
kdims=['snapshot'])
return hv.output(hmap, holomap='scrubber')
make_animation(snapshots, 5)
As can be seen in the above, the simulation stabilizes after some initial jiggling. This is due to the fact that the spring lengths set initially induce back and forth oscillations of the structure.
Random network of springs¶
Of course, this framework can be applied to different situations, such as a randomly connected network of springs.
In [9]:
n_particles = 10
x = np.c_[np.random.uniform(0.05, 0.3, size=n_particles),
np.random.uniform(0.05, 0.3, size=n_particles)]
n_springs = 20
spring_anchors = np.c_[np.random.choice(np.arange(n_particles), size=n_springs),
np.random.choice(np.arange(n_particles), size=n_springs)]
spring_stiffnesses = np.ones((n_springs,))
spring_lengths = np.random.uniform(low=0.1, high=0.2, size=n_springs)
masses = np.ones_like(n_particles)
snapshots = forward(x, steps=100, dt=0.5, damping=0.2)
In [10]:
plot_points(snapshots[-1], spring_anchors)
Out[10]:
In [11]:
make_time_plots(snapshots, 0)
Out[11]:
In [12]:
make_animation(snapshots, 5)
A structured network of springs¶
What about a structured network of springs?
In [13]:
x = np.c_[[0, 1, 2, 3, 4, 0, 1, 2, 3, 4],
[0, 0, 0, 0, 0, 1, 1, 1, 1, 1]].astype(np.float)
spring_anchors = []
for node in range(4):
# right neighbor bottom
spring_anchors.append(np.array([node, node + 1]))
# right neightbor top
spring_anchors.append(np.array([node + 5, node + 6]))
# upward diagonal right
spring_anchors.append(np.array([node, node + 6]))
# downward diagnoal right
spring_anchors.append(np.array([node+1, node + 5]))
# top down
spring_anchors.append(np.array([node, node + 5]))
spring_anchors.append(np.array([4, 9]))
spring_anchors = np.array(spring_anchors)
n_springs = spring_anchors.shape[0]
spring_stiffnesses = np.ones((n_springs,))
spring_lengths = np.ones((n_springs,))
masses = np.ones_like(x.shape[0])
snapshots = forward(x, steps=100, damping=0.1)
In [14]:
plot_points(snapshots[0], spring_anchors)
Out[14]:
In [15]:
make_time_plots(snapshots, 0)
Out[15]:
In [16]:
make_animation(snapshots, 5)
Adding some gravity and fixed nodes¶
One way we can make this system more interesting is to add some more forces to it, set some nodes to be fixed and let it move according to its dynamics.
In [17]:
def compute_gravity_forces(x, masses, g):
"""Adds a vertical force to node with mass m."""
forces = np.zeros_like(x)
forces[:, 1] -= masses * g
return forces
def forward_with_gravity_and_fixed_nodes(x, steps=200, dt=0.5, damping=0.1, g=0.1, fixed_nodes=None):
new_x = x.copy()
new_v = np.zeros_like(x)
snapshots = []
snapshots.append(new_x.copy())
for t in range(steps):
forces = compute_spring_forces(new_x, spring_anchors, spring_lengths, spring_stiffnesses)
forces += compute_gravity_forces(new_x, masses, g)
if fixed_nodes is not None:
forces[fixed_nodes] = 0
new_x, new_v = time_step_integration(new_x, new_v, forces, masses, dt=dt, damping=damping)
snapshots.append(new_x.copy())
return snapshots
In [18]:
x = np.c_[[0, 1, 2, 3, 4, 0, 1, 2, 3, 4],
[0, 0, 0, 0, 0, 1, 1, 1, 1, 1]].astype(np.float)
spring_anchors = []
for node in range(4):
# right neighbor bottom
spring_anchors.append(np.array([node, node + 1]))
# right neightbor top
spring_anchors.append(np.array([node + 5, node + 6]))
# upward diagonal right
spring_anchors.append(np.array([node, node + 6]))
# downward diagnoal right
spring_anchors.append(np.array([node+1, node + 5]))
# top down
spring_anchors.append(np.array([node, node + 5]))
spring_anchors.append(np.array([4, 9]))
spring_anchors = np.array(spring_anchors)
n_springs = spring_anchors.shape[0]
spring_stiffnesses = []
spring_lengths = []
for i in range(n_springs):
a, b = spring_anchors[i]
x_a, x_b = x[a], x[b]
dist = x_a - x_b
length = np.linalg.norm(dist)
spring_lengths.append(length)
spring_stiffnesses.append(10 * length)
masses = 2 * np.ones_like(x.shape[0])
snapshots = forward_with_gravity_and_fixed_nodes(x, steps=300, dt=0.25, fixed_nodes=[5, 0], g=0.1, damping=0.1)
plot_points(snapshots[-1], spring_anchors)
Out[18]:
In [19]:
make_time_plots(snapshots, 4)
Out[19]:
In [20]:
make_animation(snapshots, 15)
That's the end of this notebook. I hope you've enjoyed looking at these jiggly things as much as I do.
This post was entirely written using the Jupyter Notebook. Its content is BSD-licensed. You can see a static view or download this notebook with the help of nbviewer at 20191029_springsTest.ipynb.
