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)