Daniel Piker’s Kangaroo plug-in for Grasshopper allows us to set up and run various kinds of physical simulation within Grasshopper. Although we will not cover Kangaroo itself on this course, in this session we will look at a way of writing our own Kangaroo-like simulation code, allowing us to write custom simulators should the need arise and also understand how Kangaroo is (probably) working behind the scenes.
A full multi-purpose physics engine is probably a bit beyond what we are capable of producing in an hour, so instead we will focus on a specific physical phenomenon – gravitation. We will produce code that simulates a set of particles, all of which have a gravity-like attraction towards each other.
As usual, we’ll first set up a C# component with the inputs we need, which are:
pts – List Access – Point3d
dt – Item Access – double
reset – Item Access – bool
Enter into the pts input a set of points which will define the starting positions of all of our particles.
Before we begin coding, it is useful to understand a little more about Object-Orientated Programming, which is so-called because it involves the use of Objects. An object is a group of different bits of data (properties) and processes (methods) that represent a single conceptual entity. C# is an object-orientated language and we have been using objects all along – the points, lines, curves and so on that we have been creating and manipulating are all objects, as are the more abstract entities such as Lists and the random number generator we used in a previous example.
As each particle in our simulation will need to have several different properties, we will create a class in order to encapsulate all of that data into a separate object for each particle.
We will write our Particle class in the ‘custom additional code’ section of the scripting component and it will go a little something like this:
That’s a fairly big chunk of code, so we’ll go through it bit-by-bit. First of all we declare the class definition by writing:
The keyword for creating classes is, obviously ‘class’, which is followed by the name we want to give the class; in this case ‘Particle’. Everything inside the curly brackets that follow will be part of this class. The order these things are in does not particularly matter, but we will start by declaring the properties of the class:
This works similarly to declaring a local variable inside a function, however we have added the ‘public‘ keyword to the front – this means that these properties will be accessible from outside the class. If we instead wanted to make it impossible for these properties to be interfered with from outside the class itself, we could instead use the keyword ‘private‘.
We have given the Particle class three properties – ‘position’ is a Point3d and represents the current position of the particle, ‘velocity’ is a Vector3d and represents the current velocity of the particle and ‘trail’ is a list which will be used to hold all previous positions of the particle.
Having declared these three properties, we will now add functions which we can use to manipulate them. Firstly, we will create a constructor. A constructor is a function called when we create an instance of the class and whcih is used to set up the initial state of the object. Whenever we have previously use the new keyword we have been initialising an object and in the process calling some kind of constructor function. We can pass in information to the constructor to tell it how we want the object to be initialised – for example previously when we have created a new Line object and passed in two points, those were then being consumed by a constructor function which used them to set up the geometry of the line. We can have multiple constructor functions with different arguments to allow us to set up objects in different ways, however in this case we will only use the one, which takes in a Point3d:
Constructors are declared using the name of the class as a function name (in this case ‘Particle’) and omitting a return type. The inputs to the constructor are declared exactly the same way as the inputs to a function – in this case we declare an argument called ‘startPosition’ of type Point3d. Within the constructor function we set the Particle’s position property to the startPosition passed in and also initialise the velocity as a new Vector3d (which, without inputs, will start as (0,0,0)) and initialise the trail property as a new List of Point3ds. This will then set up each particle ready for use.
Note that above the constructor function is some text denoted by three forward slashes. This is a special kind of comment that, when placed just before a function definition can be used to add metadata about that function and its parameters. In any decent C# IDE, these special comments will pop up interactively whenever we use the function to give a description of how to use the function and what the input parameters represent and can also be used to automatically generate documentation (such as the RhinoCommon docs, which is entirely derived from these kind of comments in the source code). Sadly, the C# editor in Grasshopper is not a decent C# IDE, but it is still worth getting into the habit of adding these anyway.
Next, we’ll add some methods. These are functions on the object that can be called to manipulate it in various ways. We will write a function to update the position of the Particle based on its own current properties, called ‘Move’:
This function does not return anything (hence the ‘void‘ return type) and takes in a double called ‘dt’ (short for ‘delta time’). This variable represents (nominally, at least) the time that has passed since the last time the Particle was moved and is used as a scaling factor for the particle’s velocity. So if, for example, the current velocity was 1m/s and we passed in a dt of 0.5s, we multiply one by the other to find out how far the particle has moved in that time step.
Before we do that, we add the Particle’s current position to the end of the trail list (so that we have a record of all previous locations) and then modify the position by adding the velocity vector multiplied by the time step.
Finally, we add a static function to deal with the attraction between particles. Static functions do not act on a specific instance of a class and nor do they require an instance of the class in order to be called – you can access them just by typing the class name followed by a ‘.’ and the method name (i.e. ‘Particle.Attract(…)’). We have used static functions before as well – most RhinoCommon classes have some static functions, such as the function to create a curve through a set of points that we have used previously (and will also do later on). Note that execution-wise there is no particular benefit to making this function a static function on the Particle class rather than an ordinary function in our script component – it is merely a way of organising your code (and I’m doing it here mainly to demonstrate where the static functions on other object types come from).
The function itself takes in two particles and applys a pseudo-gravitational attraction between them over a period of time ‘dt’. Newton’s law of universal gravitation is:
F = G(m1*m2)/r^2
In other words; the gravitational force between two masses is equal to the gravitation constant G times the mass of the first object times the mass of the second divided by the square of the distance between them. In our case we don’t care too much about being exact, so I’ve assumed that m1, m2 and G are all equal to one (it’s our own little universe, we can do what we like) and excluded them. To find the acceleration we need to apply to each particle we then multiply that force by the timestep.
So, to begin with we find the vector that points from our first particle to our second. We then find the length of that vector (the distance between our points). Before we continue, we check to make sure that this distance is greater than 0 (otherwise, we would get a divide-by-zero error).
We then unitize the direction vector and multiply that by the timestep over the distance between the points squared – this gives us the acceleration vector to apply to each point.
[Note: With scripts like this that will run in real-time, optimising the code so that it runs quickly and responsively becomes more important. One optimisation that could be applied here (but that I haven't to keep things clear) would be to eliminate the separate call to Unitize() and instead just increase the power applied to the distance in the subsequent calculation by one (i.e. AtoB = dt/Math.Pow(distance, 3)). Since the Unitize() function simply divides a vector by its own length this will result in the same answer, but will avoid having to recalculate the Vector length as part of the Unitize() operation.]
Finally we apply the acceleration to the particle velocities. Since the vector points from particleA to particleB, we add it to particleA’s velocity to accelerate it towards B, but subtract it from particleB’s velocity to move it towards A.
That concludes all the code needed for our particle class – now let’s use it!
First, we’ll need a persistant way of storing our particles. Just above the Particle class definition, within the ‘Custom Additional Code’ section, add:
This list will be used to store all of our active particles. The beauty of declaring this outside of any particular function is that it will become a property of the script component class itself and as such will not go out of scope once our RunScript function has ended. This means that the data will be retained – the next time we trigger an update on the component and RunScript is called, the particles list will still contain exactly the same thing it contained the last time RunScript finished. This allows us to have a set of Particles that persist between firings of the component. Each time the component is triggered, we will update the particles in this list slightly.
Now, finally, on to the main RunScript function, which will look like this:
To begin with, we need a way of setting up the initial state of the simulation. Whether this bit of code will be run or not depends on whether the ‘reset’ input is set to true. If so, we will use this to restore the simulation to it’s starting condition – otherwise we will allow it to progress using its current state.We first overwrite the particles list with a new, empty, list of particles. For each point that we’ve input we create a new particle by using the new keyword. We pass in the point as an argument to the constructor that we wrote earlier, which initialises the Particle position to be at that point. We then add the particle to the particles list.
When the reset option is false, we instead run the actual simulation.
The first thing we need to simulate is the motion of the particles. To do this we simply iterate through all of the particles and call the Move function, passing in the specified timeStep. This will update the position of each particle based on its own individual velocity.
The next part is slighty more complicated. We want to adjust the velocity of each particle to attract it towards all other particles. To do this we need to call the Attract function once for each possible pairing of particles. To ensure that each pairing only gets processed once we use a pair of nested loops. The outer loop extracts each particle in the particles list apart from the very last one one-by-one and uses that as particleA. The inner loop extracts each particle in the list after the position of particleA and uses that as particleB. By doing this we ensure that each particle pairing is only used once – we do not need to compare particleA to anything before it in the list because those pairings will have already been covered.
That’s it – our simulation is now ready to run. However, we won’t yet be able to see it running, so the final bit of code we need to write outputs the current position of each particle and also uses the trail stored on the particle to create an interpolated curve which describes it’s past path. We then output these to the A and B outputs of our component.
Ta da! Our code is complete! Now close the component and you should see that… nothing has happened. This is because, in order for our simulation to work it needs to be run over and over and over again. However, grasshopper will only run the code by default when one of the inputs changes.
To get the component to constantly recalculate, we need to use a Timer component (Params/Util). Add one of these to the canvas. Right-click on it and set the interval to 20ms. Then drag the output from the right hand side of the timer and connect it to the C# script component (the Timer works differently to most other components – you do not have to connect it to a specific input, it works on the component itself).
The timer component will now cause our C# component to automatically re-fire every 20 milliseconds. Toggle reset to true, then to false and you should see the simulation play out before your eyes, with the particles tracing curling paths around one another through space.
- You can use this same basic approach to simulate pretty much any physical process you like, simply by modifying the interaction between the particles. If you wanted to model a bunch of objects all bouncing off of one another, you could replace the Attract() function with code to detect collisions between particles and calculate their resultant velocities based on conservation of momentum. If you wanted to study boid-like flocking behaviour, replace it instead with code that steers each particle with respect to those around it. The possibilities are endless!
- You could simulate the same behaviour without the timer and the mechanisms used to persist data simply by looping inside the component instead. This will be faster (because Rhino will not be constantly re-drawing the output) but has the disadvantages that you will not be able to see or interact with the process while it is running. It is therefore what we structural engineers call ‘less fun’.
- Experiment with tweaking the input parameters. You will probably see that adjusting the timeStep, as well as speeding up or slowing down the motion of the particles, also has an effect on the final result you get, even with the same set of starting locations. This is due to a simplification in the way we are modelling the motion of the particles. We are treating each particle as being accelerated in little instantaneous ‘jumps’, moving in straight lines in between. In reality, the gravitational attraction is being applied constantly, leading to a slight inaccuracy in the positions of the particles at each step. Because the subsequent acceleration is dependant on the distance between particles, this leads to further little inaccuracies which over time add up to quite major errors. To calculate the motion more accuractely is a slightly more complicated process, although as designers rather than astronomers we may decide we can live with the simpler version.
- I have used a similar orbital simulation when working on the design of the ArcelorMittal Orbit, with the path of a single particle moving around fixed attractor points being used to determine the central ‘spine’ around which the diagrid geometry was generated. We discovered that the simple form of this was a little too sensitive to initial conditions (a minor tweak to one of the control points would completely change the entire path) and so the model was refined to allow the particle to ‘steer’ itself towards guidance points which allowed us to better adjust the geometry to meet artistic and structural drivers.
The example file for this session can be downloaded here: