Ropes

This example simulates ropes as nodes connected by slightly deformable sticks.

# %pip install --quiet phiflow
from phi.jax.flow import *
# from phi.flow import *  # If JAX is not installed. You can use phi.torch or phi.tf as well.
from tqdm.notebook import trange

We start by scattering nodes as a jittered grid and setting some of them as fixed points while the rest will be able to move.

grid = CenteredGrid(0, x=20, y=20, bounds=Box(x=1, y=1))
x = pack_dims(grid.points, 'x,y', instance('nodes'))
x += math.random_uniform(x.shape) * .01
fixed_indices = vec(nodes=[399, 19, 199, 239, 0])
fixed = math.scatter(expand(False, instance(x)), fixed_indices, True)
plot(PointCloud(x, fixed), size=(4, 3));

C:\PhD\phiflow2\examples\particles\../..\phi\vis\_matplotlib\_matplotlib_plots.py:167: UserWarning: This figure includes Axes that are not compatible with tight_layout, so results might be incorrect.
  plt.tight_layout()  # because subplot titles can be added after figure creation

We represent our rope system using a Graph whose edges encode the stick lengths. The sticks are represented by the graph's edge matrix, which we store in sparse COO format here.

deltas = math.pairwise_differences(x, max_distance=grid.dx.mean * 1.1, format='coo')
distances = math.vec_length(deltas)
graph = geom.graph(x, distances, {})
plot(graph);

C:\Users\phili\AppData\Local\Temp\ipykernel_20812\2855547231.py:2: DeprecationWarning: phiml.math.length is deprecated in favor of phiml.math.norm
  distances = math.vec_length(deltas)

Now let's define the simulation! After handling the gravity and velocity updates, we perform a fixed number of relaxation steps.

Each relaxation step shifts the node positions in order to reduce stretching or squishing of the sticks. To do this, we first compute the stick centers and directions. For an unstretched stick, the node position must be half the stick's length away from its center along the stick's direction. For nodes with multiple connected sticks, we compute the ideal node position for each stick, and then average the so-obtained positions.

@jit_compile
def step(graph: geom.Graph, v, dt=1., gravity=vec(x=0, y=-0.01), relaxation_steps=50):
    v += gravity * dt
    x = graph.center + math.where(fixed, 0, dt * v)
    for _ in range(relaxation_steps):
        deltas = math.pairwise_differences(x, format=graph.edges)
        stick_centers = x + .5 * deltas
        dist = math.norm(deltas)
        stick_directions = deltas / (dist + 1e-5)
        next_x = stick_centers - stick_directions * .5 * graph.edges
        next_x = math.mean(next_x, dual)
        x = math.where(fixed, x, next_x)
    v = (x - graph.center) / dt
    return geom.graph(x, graph.edges), v

v0 = math.zeros_like(x)
graph_trj, v_trj = iterate(step, batch(time=50), graph, v0, substeps=2, range=trange)

  0%|          | 0/50 [00:00<?, ?it/s]

show(graph_trj, animate='time')

<Figure size 640x480 with 0 Axes>