#### Basic physics for simulating fabrics

A common method to simulate cloth for computer graphics and the movie industry is to use a so called spring and dampener system. This method is based on simulating the fabrics by splitting them up on a grid and simulate the movement of the cloth at each point (intersection) on the grid. By enforcing forces and constraints that act on these points you can simulate many different types of fabrics, and by using a finer resolution grid you get more natural and realistic simulations.

By varying different constants during the simulation you can simulate many different types of fabrics. In the example below we are going for a fairly stiff and smooth material.

Each point on the grid can be simulated using Newtonian physics to update the velocity and position of the points. However, in order for the grid to behave as a fabric and not just a collection of points we need to add constraints and forces that act on these points. Typically this is done with the metaphor of springs and dampeners. We imagine that there exists a spring-like force between nearby points that applies a force to keep them at a set distance. Likewise we imagine that there exists dampeners in additions to the springs that simulates friction, getting rid of the unwanted oscillations that one would otherwise get from a perfect spring.

Consider the grid with springs and dampeners in the illustration above. We can here define the force acting on each point as a function of the extensions of the springs and the relative velocities of the two points. If we store the position and velocities that act on the grid points in the arrays P and V we can compute the forces that act on the points as follows. For the case of a single connection between grid cell [i,j] and [i+1,j]:

where k_s and k_d are the spring constant and dampening constant, respectively and l is the spring length at rest.

In order to efficiently calculate the forces over the whole grid, we use array operations to perform the calculations above in a single step. In numpy this can be formulated as follows for the spring part of the calculation:

# skips the last point

pos1 = pos[:-1, j, :]

# skips the first point

pos2 = pos[1:, j, :]

# vector from pos1 to pos2

v12 = pos2 - pos1

# length of vector

r12 = np.linalg.norm(v12, axis=1)

f = (v12.T * (r12 - size)).T * ks

self.force[:-1, j, :] += f

self.force[1:, j, :] -= f

Note how the code above computes the force for all ”horizontal” connections, and updates the forces on both points that are connected together.

To compute the updated positions and velocities we can use any of the standard physics integrators such as an Euler integrator (worst), Verlet integrator (better) or even a Runge-Kutta integrator (RK4, best). In our case we pick a simple Euler integrator since it uses a formulation that is easier to translate to Runge-Kutta RK4 in the end.

#### Structural forces

If we only consider connections between the four direct neighbours (left, right, up, down) then we will not get a very life like material. What we need to do is to preserve forces in the plane of the object (structural forces), and shearing forces as well as to resist bending forces.

The structural forces can be given directly with the four connections to the direct neighbours. These connections are illustrated with the black dashed lines below.

For the shearing forces, we need to add constraints along the diagonal, again by adding springs and optionally dampeners. This is illustrated by the red arrows below.

Finally, in order to prevent the fabric from bending too much we add springs and optionally dampeners along connections of length 2. This alone allows the springs and dampeners to preserve rotational forces by adding a force only whenever the distance between two points are not exactly twice the resting length of the springs, this happens if and only if the grid is bent. Thus creating a force that straightens the grid. These connections are illustrated by the green arrows below.

Note that the resting length of each spring is calculated as per the length of the corresponding connection in the the original grid, ie. the diagonals have sqrt(2) times the length of the direct connections. By changing the resting length (and connectivity) of the connections we can alter the default shape of the object. By changing the spring constants and dampening constants we can alter the type of material that is simulated.

#### Accelerating the simulation using TensorFlow

To speed up the calculations we will use TensorFlow to perform all array calculations. To get started with tensorflow take a look at one of the many tutorials that are available. In the code below we use the original tensorflow method of first building up a graph that describes the calculations to perform, as opposed to using the slower tensorflow eager mode that perform the calculations on the go.

We start by defining the position and velocities of the grid as tensorflow variables.

args = {"trainable": False, "dtype": tf.float32}

pos_t = tf.Variable(pos, name="pos", **args)

vel_t = tf.Variable(vel, name="vel", **args)

where pos and vel are pre-existing numpy arrays that contain the starting positions and starting velocity.

We define useful constants for the calculations such as step time, gravity etc. as tensor constants:

mass = tf.constant(1.0, name="mass")

dt = tf.constant(2e-3, dtype=tf.float32, name="dt")

gravity = tf.constant(np.array([0.0,-9.81,0.0]), dtype=tf.float32, name="g")

size = tf.constant(0.1, dtype=tf.float32, name="size")

...

We can calculate the internal forces that act on the cloth by converting the elementwise operation from before into tensor operations. We accumulate all the different force calculations into a list of tensors, and perform a final summation of them as a last step.

forces = []

# Spring forces for i+1

pos1 = pos[:-1, :, :]

pos2 = pos[1:, :, :]

v12 = pos2 – pos1

r12 = tf.norm(v12, axis=2)

f = (v12 * tf.expand_dims((r12 – size), axis=2)) * ks

f_before = -tf.pad(f, tf.constant([[1,0],[0,0],[0,0]]))

f_after = tf.pad(f, tf.constant([[0,1],[0,0],[0,0]]))

forces.append(f_before)

forces.append(f_after)

…

# Dampening for i+1

vel1 = vel[:-1, :, :]

vel2 = vel[1:, :, :]

f = (vel2 – vel1) * kd

f_before = -tf.pad(f, tf.constant([[1,0],[0,0],[0,0]]))

f_after = tf.pad(f, tf.constant([[0,1],[0,0],[0,0]]))

forces.append(f_before)

forces.append(f_after)

…

total_force = tf.add_n(forces)

Similarly to the example above we add the calculations for:

- spring forces for neighbours on the same row/column: i+1, i+2, j+1, j+2,
- spring forces for neighbours on diagonals: (i+1, j+1), and (i-1, j+1)
- dampening forces for direct neighbours: i+1, j+1

We can also add collision forces with a ball and the ground in order to make for a more interesting simulation:

ball_center = tf.constant(np.array([0,0,0]), dtype=tf.float32, name="ball")

ball_radius = tf.constant(1.0, dtype=tf.float32, name="rad")

V = pos - ball_center

r = tf.norm(V, axis=2)

p = tf.maximum(tf.constant(0, dtype=tf.float32, name="0"), ball_radius - r)

r = tf.reshape(r, r.shape.dims+[tf.Dimension(1)])

p = tf.reshape(p, p.shape.dims+[tf.Dimension(1)])

force = V / r * p * tf.constant(1e5, dtype=tf.float32)

Note how we can reformulate the problem of collision detection into computing how far into the ball a point is and apply an outwards force, using a MAX operation to ensure that points that are not inside the object (ball_radius – r is negative) are not affected. This allows us to do collision detection without any conditional operators which would have been slow otherwise.

Finally, we add collisions with the ground in a similar manner and add gravity that affects all points.

The update step is done as a naive Euler implementation (for now).

self.delta_vel = self.force4 * self.dt / mass

self.delta_pos = self.vel_t * self.dt

self.step = tf.group(

self.vel_t.assign_add(self.delta_vel),

self.pos_t.assign_add(self.delta_pos))

When generating the animation below we alternate between a number of physics steps and extracting the position data from GPU memory to CPU memory and visualising the mesh:

for i in range(100): session.run(self.step)

self.pos = self.pos_t.eval(session=session)

self.draw()

If we where to use a more advanced integrator such as RK4 we could have a larger step-size without introducing instabilities to the simulation. That would overall serve to speed up the simulation and allow us to get away with fewer calls to TensorFlow.

#### Execution time

To see how effective our tensorflow implementation was we measured the time per single update step, ie. calculating the forces and updating the velocity and positions one time. We measure this time as a function of the total array size (width * height) since we know that the total execution time scales by the total number of points or the square of the grid-resolution.

As we can see in the graph below we have a fairly large constant offset in time that makes the number of points used barely have any effect below 40.000 points (200 in width and height). This is most likely caused by the time needed to start the GPU kernels containing the calculations above, which puts a limit to how useful this method is when simulating many smaller objects.