Differential Growth

Differential growth is a process where different parts of a structure grow at different rates, leading to complex forms. After watching this, I decided to try it in Grasshopper. Here, differential growth mimics this behavior by applying rules such as Repulsion to avoid crowding or overlapping, Cohesion to keep parts connected or within a range, and Insertion to add new elements when a part stretches too far. The algorithm I present here simulates a differential growth process starting from a set of points along a curve. All agents apply repulsion forces to each other if they are within a certain distance. Each agent checks its neighbors: it pushes or pulls them. Additionally, when two neighbors become too far apart, the code inserts a new (point) agent. Over time, this interaction leads to growing patterns shown below:

I developed this Grasshopper Python code in Rhino 8 (Python 3). This means it will NOT work in Rhino 7. Crv is the initial curve for the starting agents. I is the distance threshold above which a new agent is inserted. R is the range within which global repulsion is applied. K is the scaling factor for force accumulation. Q is the damping coefficient (controls how fast they slow down). N is the maximum number of agents. When Run is triggered, it triggers one iteration of the simulation. I attached a Grasshopper Timer component to repeat it. The outputs are the points and the vectors. This implementation makes some optimizations, but it still checks every boid with every other. So it becomes very (really) slow when you increase the N value too much. Below is the Python code inside the Grasshopper definition. You can find the definition at the bottom of this page.

import rhinoscriptsyntax as rs
import math
insertDistance = I
repulsionDistance = R
k = K
q = Q
def vector_add(a, b):
	return [a[0]+b[0], a[1]+b[1], a[2]+b[2]]
def vector_subtract(a, b):
	return [a[0]-b[0], a[1]-b[1], a[2]-b[2]]
def vector_scale(v, s):
	return [v[0]*s, v[1]*s, v[2]*s]
def vector_length(v):
	return math.sqrt(v[0]**2 + v[1]**2 + v[2]**2)
def vector_unitize(v):
	length = vector_length(v)
	if length == 0:
		return [0,0,0]
	return [v[0]/length, v[1]/length, v[2]/length]
def vector_distance_sq(a, b):
	return (a[0]-b[0])**2 + (a[1]-b[1])**2 + (a[2]-b[2])**2
class Boid():
	def __init__(self, pos, vel):
		self.pos = list(pos)
		self.vel = list(vel)	
	def updateVel(self, i):
		forces = [0,0,0]
		for boid in flock:
			if boid != self:
				dis = vector_distance_sq(boid.pos, self.pos)
				if dis < repulsionDistance ** 2:
					dif = vector_subtract(boid.pos, self.pos)
					dif = vector_unitize(dif)
					dif = vector_scale(dif, -1 / dis)
					forces = vector_add(forces, dif)
		n_boid = flock[i-1]
		p_boid = flock[(i+1) % len(flock)]
		p_dis = vector_distance_sq(p_boid.pos, self.pos)
		p_dif = vector_subtract(p_boid.pos, self.pos)
		p_dif = vector_unitize(p_dif)
		p_dif = vector_scale(p_dif, 1 / p_dis)
		if p_dis < (insertDistance * 0.5) **2:
			p_dif = vector_scale(p_dif, -1)
		forces = vector_add(forces, p_dif)
		n_dis = vector_distance_sq(n_boid.pos, self.pos)
		n_dif = vector_subtract(n_boid.pos, self.pos)
		n_dif = vector_unitize(n_dif)
		n_dif = vector_scale(n_dif, 1 / n_dis)
		if n_dis < (insertDistance * 0.5) **2:
			n_dif = vector_scale(n_dif, -1)
		forces = vector_add(forces, n_dif)
		forces = vector_scale(forces, k)
		self.vel = vector_add(self.vel, forces)
	def updatePos(self):
		self.pos = vector_add(self.pos, self.vel)
		self.vel = vector_scale(self.vel, q)
def check_insert_particles(flock):
	i = 0
	while i < len(flock):
		boidA = flock[i]
		boidB = flock[(i + 1) % len(flock)]
		dis = vector_distance_sq(boidA.pos, boidB.pos)
		if dis > insertDistance**2:
			pos = vector_scale(vector_add(boidA.pos, boidB.pos), 0.5)
			vel = vector_scale(vector_add(boidA.vel, boidB.vel), 0.5)
			flock.insert(i+1, Boid(pos, vel))
		else:
			i += 1
def reset_flock():
	points = rs.CurvePoints(Crv)
	points = rs.CullDuplicatePoints(points)
	center = rs.CurveAreaCentroid(Crv)[0]
	for point in points:
		vel = vector_subtract(point, center)
		vel = vector_unitize(vel)
		boid = Boid(point, vel)
		flock.append(boid)
if "ready" not in globals() or Reset:
	ready = True
	flock = []
	reset_flock()
if ready and Run:
	if len(flock) < N:	
		check_insert_particles(flock)
	for i, boid in enumerate(flock):
		boid.updateVel(i)
	for boid in flock:
		boid.updatePos()
pos = [rs.CreatePoint(*boid.pos) for boid in flock]
vel = [rs.CreateVector(*boid.vel) for boid in flock]