Verlet Integration

Physics Simulation in Clojure

Posted by: Maarten Metz on 2018-03-04
ClojureFunctional

Introduction

In this post we’re going to dive into physics simulation with Clojure. The ideas in this post were mostly inspired by this paper by Thomas Jakobsen.

My motivations for working on this toy project were:

  • Experimenting with a very mutable domain in a language with immutability at its core
  • Having some fun by making a lot happen on screen without writing lots of code
  • Coding using emacs & CIDER (The Clojure Interactive Development Environment that Rocks)

The source code can be found here. The result:

One main- to bind them all

This code uses the Quil library to render points and lines. In the main- function, the sketch macro from the Quil library is used to bind the physics domain to the UI.

(defn -main []
  (quil/sketch
    :host           -main
    :title          "Verlet Integration"
    :size           [width height]
    :setup          setup
    :update         update-state
    :draw           draw
    :key-pressed    key-pressed
    :mouse-pressed  mouse-pressed
    :mouse-dragged  mouse-dragged
    :mouse-released mouse-released
    :features       [:exit-on-close]
    :middleware     [quil-mw/fun-mode]))

The sketch macro basically asks the developer to configure which ‘handler functions’ it should call when certain events occur. For instance:

  • a key is pressed -> call key-pressed
  • the mouse is pressed -> call mouse-pressed
  • the state of your program has to be updated -> call update-state
  • it’s time to draw the new state -> call draw
  • etc.

Because I configured Quil to run in functional mode (quil-mw/fun-mode), Quil:

  • Uses the return value from the :setup function as the initial state
  • Passes the state to each handler function to update it
  • Additionally passes a keyboard- and mouse event argument to keyboard - and mouse handler functions.

The handler function signatures therefore look like this:

(defn draw          [state]       ...)
(defn update-state  [state]       ...)
(defn key-pressed   [state event] ...)
(defn mouse-pressed [state event] ...)

World state

What does a snapshot of this world state look like? key-pressed gives some hints:

(defn key-pressed
  [state event]
  (let [raw-key   (:raw-key event)
        new-state (cond
                   (= \b raw-key) (load-world (load-a-file "blocks.edn"))
                   (= \c raw-key) (load-world (load-a-file "cloth.edn"))
                   (= \p raw-key) (load-world (load-a-file "points.edn"))
                   (= \s raw-key) (load-world (load-a-file "sticks.edn"))
                   (= \i raw-key) (show-info-message)
                   (= \q raw-key) (quil/exit)
                   :else state)]
    new-state))

So on pressing a key, a certain world is loaded. Let’s look into the ‘sticks’ world:

{:points {:p0 {:x  10 :y   1 :oldx   0 :oldy   0}
          :p1 {:x 100 :y 100 :oldx 100 :oldy 100}
          :p2 {:x   0 :y 500 :oldx   0 :oldy 525}
          :p3 {:x  25 :y 475 :oldx   0 :oldy 525}
          :p4 {:x 250 :y 200 :oldx 250 :oldy 200 :pinned true}
          :p5 {:x 350 :y 100 :oldx 350 :oldy 100}}
 :sticks [{:links  [:p0 :p1] :length 138}
          {:links  [:p2 :p3] :length  40}
          {:links  [:p4 :p5] :length 140}]}

Nothing special here. Just a map containing :points and :sticks. We’ll dive into the particulars of points and sticks soon.

Physics Simulation

This physics simulation deals with Points, and Sticks connecting Points. The Points are the things that seem to have direction and speed, are influenced by gravity and lose speed because of friction or because of bouncing against world borders. The Sticks try to keep their 2 Points apart according to the configured Stick length. The simulation loop boils down to:

  • Update Points
  • Apply Stick constraints
  • Apply world constraints

In code:

(defn update-state
  [state]
  (->> state
       (update-points)
       (apply-stick-constraints)
       (apply-world-constraints)))

This code can be read like this: ‘with state, first update-points, then apply-stick-constraints, and finally apply-world-constrains’.

Points

Points are the main abstraction in this code. I decided to use a record to name them:

(defrecord Point [x y oldx oldy pinned])

A Point stores its current coordinates and the coordinates it had in the previous world state. On top of that it has a pinned property which indicates if a Point is pinned in space and - as a result - stays on the same coordinate. A pinned Point can be unpinned by clicking it with the mouse pointer.

The update-point function calculates the velocity of the Point and adds some gravity in the mix:

(defn update-point
  [{:keys [x y oldx oldy pinned] :as point}]
  (if pinned
    point
    (let [vx (velocity x oldx)
          vy (velocity y oldy)]
      (->Point (+ x vx) (+ y vy gravity) x y pinned))))

Reading superficially, this code says: if the Point is pinned then return the same Point. Otherwise, calculate it’s velocity (based on the current - and previous x & y coordinates), add some gravity in the y direction and return the newly calculated Point.

The function uses destructuring to name all the arguments of the incoming Point. With destructuring, you can bind the values in a data structure without explicitly querying the data structure. So instead of getting each and every value out of Point and binding it to a ‘new’ name…

(let [x      (:x      point)
      y      (:y      point)
      oldx   (:oldx   point)
      oldy   (:oldy   point)
      pinned (:pinned point)]
  ...)

… you can very effectively tell the function that it will receive a map {} with keys x y oldx oldy pinned and that this function should have those variables available under the same name:

[{:keys [x y oldx oldy pinned] :as point}]

The code in my repository has some additional type hints, because I wanted the simulation to run smooth on my laptop and wanted to learn a bit more about type hinting and compiler optimizations.

And that’s it with regard to Points. In a couple of lines of code Points are already moving and reacting to gravity in the simulation.

World constraints

It’s time for Points to meet the harsh reality of life. Walls are harder than Points and therefore Points should bounce off of them:

(defn hit-floor?       [y] (> y height))
(defn hit-ceiling?     [y] (< y 0))
(defn hit-left-wall?   [x] (< x 0))
(defn hit-right-wall?  [x] (> x width))


(defn apply-world-constraint
  [{:keys [x y oldx oldy pinned] :as point}]
  (let [vxb (* (velocity x oldx) bounce)
        vyb (* (velocity y oldy) bounce)]
    (cond
      (hit-ceiling?    y) (let [miry (- oldy)]
                             (->Point (+ oldx vxb) (+ miry (- vyb)) oldx miry pinned))
      (hit-floor?      y) (let [miry (+ height height (- oldy))]
                             (->Point  (+ oldx vxb) (- miry vyb) oldx miry pinned))
      (hit-left-wall?  x) (let [mirx (- oldx)]
                             (->Point (+ mirx (- vxb)) (+ oldy vyb) mirx oldy pinned))
      (hit-right-wall? x) (let [mirx (+ width width (- oldx))]
                             (->Point (- mirx vxb) (+ oldy vyb) mirx oldy pinned))
      ;; else: free movement
      :else point)))

To explain what happens here, a picture might help.

apply world constraints diagram

Imagine a Point moved from A to B in the last Point update. This means a Point record is persisted with its x and y being the coordinates of B and oldx and oldy being the coordinates of A. If a wall line crossed the imaginary line A-B, the Point history should be rewritten. In essence, line A-B is mirrored in the wall it hits, giving rise to another imaginary line C-D where C mirrors A and D mirrors B. The miry (mirror-y) and mirx (mirror-x) variables in the code contain C coordinates mirroring A coordinates.

The simulation will take a little velocity loss into account when bouncing at walls. Therefore, D’ is calculated by using the x and y velocities (vx and vy) multiplied by a bounce factor (leading to the vxb and vyb variables in the code).

When the world state is drawn, the Point will fly off in exactly the right direction. Additionally, because of the history rewrite, subsequent Point updates will keep moving the Point in the right direction.

The picture shows the situation when hitting the ceiling. Hitting the floor and the walls works similar. And that’s all the math and code you need to simulate Points moving in a bounded space. Now let’s add Sticks to the simulation.

Stick constraints

Not only walls are restricting Points from free movement; Sticks also constrain them. A Stick connects 2 Points and has a configured length. The goal of the Stick constraint is to move the Points at the end of the Stick to the configured length of the Stick.

apply stick constraints diagram

Instead of trying to calculate the solution that satisfies all Stick constraints at once, this code simply looks at one Stick at a time and ‘solves’ the constraint problem by repeatedly solving Stick constraints in isolation.

(defn apply-stick-constraint
  [{length :length :as stick}
   {p0x :x p0y :y oldp0x :oldx oldp0y :oldy pinp0 :pinned :as p0}
   {p1x :x p1y :y oldp1x :oldx oldp1y :oldy pinp1 :pinned :as p1}]
  (let [{:keys [dx dy distance]} (distance-map p0 p1)
        difference (- length distance)
        fraction   (/ (/ difference distance) 2)
        offsetX    (* dx fraction)
        offsetY    (* dy fraction)
        p0-new     (->Point (- p0x offsetX) (- p0y offsetY) oldp0x oldp0y pinp0)
        p1-new     (->Point (+ p1x offsetX) (+ p1y offsetY) oldp1x oldp1y pinp1)]
    [(if pinp0 p0 p0-new) (if pinp1 p1 p1-new)]))

apply-stick-constraint takes a Stick and 2 Points (P0 and P1) as arguments. These arguments are heavily destructured. To understand what happens precisely, the next picture might be helpful:

stick math diagram

Let P0 and P1 be two Points with a Stick S between them. P0 and P1 have already been updated by the update-point function and their new (x,y) coordinates are (40,50) and (70,90) respectively. This means:

  • dx = 30
  • dy = 40
  • distance = 50

Since the Stick between P0 and P1 was configured to have a length of 150, and we want P0 and P1 to be adjusted by the same amount, P0 should be located at Point S0 (10,10) and P1 should be located at Point S1 (100,130). This is exactly what apply-stick-constraint achieves:

  • The difference in the configured length of the Stick and the distance of P0 and P1 is 100
  • The fraction of the distance between P0 and P1 that needs to be added to P0 AND to P1 is 1. So by extending the distance P0 and P1 with exactly 50 (distance times fraction) in both directions, we’ll solve the Stick constraint.
  • This is done by using offsetX on P0x and P1x and using offsetY on P0y and P1y

Conceptually a scaled helper triangle is used to move P0 and an identical scaled helper triangle is used to move P1. The scale of this triangle is determined by the fraction which in this case is 1. These triangles are depicted in red.

Each Stick in the world is updated before a new world frame is drawn. It is very well possible that the neatly placed P0 and P1 are moved to different coordinates by another Stick constraint before the complete update is over. This sometimes results in wiggly behavior, although in most cases, the simulation is perceived to behave naturally.

Drawing

Finally, the Points and Sticks have to be drawn on screen:

(defn draw
  [state]
  (quil/background 255)
  (if (:info-message state)
    (quil/text info-message 20 20)
    (do
      (doseq [point (vals (:points state))]
        (quil/ellipse (:x point) (:y point) 7 7))
      (doseq [stick (:sticks state)]
        (let [points (:points state)
              p0     ((first (:links stick)) points)
              p1     ((last  (:links stick)) points)]
          (quil/line (:x p0) (:y p0) (:x p1) (:y p1)))))))

Start by clearing the screen. If there is an :info-message entry in the state, draw the info-message. Otherwise draw all Points in state as a quil/ellipse and draw all Sticks as a quil/line. And that’s all folks!

Conclusion

As always, it was a pleasure working with Clojure. It turns out to be easy to work in a very mutable domain with mostly pure functions, also thanks to the excellent Quil library.

My current version of cloc tells me that core.clj is 230 lines long, but it counts the function documentation as code. Therefore, I ran a new cloc command without the doc strings. It then counts 159 lines of clojure, approximately half of which are dedicated to UI interaction and the other half to physics simulation. Not bad at all!

I’m glad I decided to ‘bite the bullet’ and learn Emacs and CIDER. I now use Emacs for all my Clojure(script) development. In addition, I use Emacs for most of my text editing these days including writing and planning.

I’d like to thank Michiel Borkent a.k.a. @borkdude for reviewing an earlier version of the code and giving me very helpful feedback. Faults and not-so-idiomatic Clojure code remaining are my own.

Please share your comments, suggestions and thoughts about this blog post on twitter.com/mmz_. Thanks for reading and Happy Coding!