Hybrid programming 2

Introducing python - N-bodies

We continue to refactor another previous example, in this case, from the Multibody Dynamics session [>] to demonstrate a hybrid visual-textual programming approach.

Sun and Earth revisited

For the refactored Sun and Earth grasshopper program shown above, we reuse PyInitialise and PySave from before as well as the original setup (red) and display (yellow) blocks. However, PyPhysics needs to be reimplemented as we must now capture the notion of n-ary gravitational forces between multiple bodies.

3 Python components are created in place of PyPhysics. The first, PyGravityForce, appears to have an unconventional setup. It has 5 input parameters. pointA and pointB receive the same input positions data, except that the latter is grafted. A similar situation occurs with massA and massB.

Each input parameter is set to default ‘item access’. The python component works like any other Grasshopper component and attempts to automatically match input data. Given an input list with 2 items and an input datatree with 2 branches each storing one item (or a 2-celled column and a 2-celled row using the spreadsheet analogy), the component loops implicitly and outputs a tree with 2 branches, each storing two items (or a 2×2 table).

from Rhino.Geometry import Vector3d

vecBA = pointA - pointB
squared_dist = vecBA * vecBA

if squared_dist == 0:
    a = Vector3d.Zero
else:
    vecBA.Unitize()
    a = (G * massA * massB / squared_dist) * vecBA

The above code computes the gravity force acting on a particle based on Newton’s law $F = GmM/d^2$. A vector going from B to A is created by point subtraction and the dot product of the vector by itself gives its squared length $d^2$. Due to the automatic data matching behaviour of the component, we get the cartesian product i.e. all pairwise combinations of particles from which we can calculate vectors and mass products. vecBA is unitized and scaled by $F = GmM/d^2$ to get the gravity force vector. The code also checks for particle-self pairs and outputs a zero vector in such cases.

drag_force = K * -velocity
a = drag_force + grav_force

All gravity force vectors acting on a particle are added up and then passed to PySumForces which computes drag force and sums it with grav force. These forces are subsequently passed to PyIntegrate to compute future states of the particle as well as apply constraints. We now extend the program to be able to simulate multiple randomly positioned particles within a bounded universe.

47KB

2_Sun_Earth_Py.gh

N-body setup revisited

if enable:
    # Change vector component direction
    if position_new.X < -limit or position_new.X > limit:
        velocity_new.X *= -1
    if position_new.Y < -limit or position_new.Y > limit:
        velocity_new.Y *= -1
    if position_new.Z < -limit or position_new.Z > limit:
        velocity_new.Z *= -1

The constraints related code in PyIntegrate is modified. If a particle’s X, Y or Z coordinates falls outside a range [-limit, limit], we change the sign of its velocity’s X, Y or Z component. Here the operator *= multiplies the right operand (-1) with the left operand (velocity component) and then assigns the result to the left operand velocity_new.X. This is an example of an augmented assignment statement.

26KB

3_Nbody_Py.gh

Iteration

We now consider refactoring the n-bodies program further and reduce physics related components to a single one—PyPhysics. To do so, we no longer rely on the convenience of grafting to get cartesian product, but instead must explicitly find every pairwise combination. Here we introduce the concept of iteration—repeated execution of a block of code—for this task. For PyPhysics, positions , velocities and masses are set to list access instead of item access by default. After doing so, we iterate through these lists rather than allow the component to create its own implicit loops.

from Rhino.Geometry import Vector3d
# Create empty list to store forces, velocities and positions
forces_new = []
velocities_new = []
positions_new = []

num_particles = len(positions)
# Iterate through particles
for i in range(num_particles):
    #1) Calculate force
    point0 = positions[i]
    mass0 = masses[i]
    grav_force = Vector3d.Zero
    for j in range(num_particles):
        if i != j:
            point1 = positions[j]
            mass1 = masses[j]
            vec01 = point1 - point0
            f = G * mass0 * mass1 / (vec01 * vec01) # F = G.m.M/r*r
            vec01.Unitize()
            grav_force += (vec01 * f)
    drag_force = K * -velocities[i]
    sum_force = grav_force + drag_force
    
    # 2) Integrate
    accel = sum_force/mass0
    dv = accel * dt
    velocity_new = velocities[i] + dv
    dp = velocity_new * dt
    position_new = point0 + dp

    # 3) Constraints
    if enable:
        # Change vector component direction
        if position_new.X < -limit or position_new.X > limit:
            velocity_new.X *= -1
        if position_new.Y < -limit or position_new.Y > limit:
            velocity_new.Y *= -1
        if position_new.Z < -limit or position_new.Z > limit:
            velocity_new.Z *= -1

    forces_new.append(sum_force)
    positions_new.append(position_new)
    velocities_new.append(velocity_new)

a = positions_new
b = velocities_new
c = masses
d = forces_new

The entire code to analyse forces, integrate and apply constraints is shown above. At this stage, having not introduced Python formally in the course, we do not expect you to fully understand every statement. Nonetheless we encourage the reader, especially those with prior programming experience, to take a stab at comprehending the code. Take note of the use of the for loop here, a compound statement that takes the general form

For <variable> in <sequence>:
	<block of code>

We have an outer for loop and another nested inside it, creating loops within a loop. This is the process by which we compute all particle pairs and mass product pairs.

Practice

The itertools module contains functions for iterating over data structures. Try importing it and using the product() function

Conclusion

In this session, we presented a hybrid approach of mixing coding with visual programming. The resulting Grasshopper program is smaller (containing fewer graphic tokens) as symbolic expressions allow for more compact representations of algorithmic logic. We also avoided having clusters nested within clusters, which arguably introduces a level of opacity to the program. Mixing in python coding thus allows Grasshopper programs to scale in complexity while maintaining ease of use and readability.

At the same time, we took this opportunity to introduce Python briefly, to build some familiarity with the language before it is introduced more extensively and in greater depth later in this course. As with any language, there is an initial barrier to overcome in terms of learning its syntax and semantics. Yet it is beneficial to be exposed to different programming languages, learn their affordances, and note how they support or hinder various modes of thinking.