7. Dynamic Simulation of Spring-Mass Systems¶

In the kinematic simulations described in the previous chapters, there was no interaction between the particles. In this chapter, we will learn how to program more complex interactions between objects through dynamic simulations.

This may sound difficult, but by using the concepts we have learned so far, such as class inheritance and object composition, we will see that it is possible to build complex systems with simple combinations of elements.

Left click: create a particle
Right click: create a fixed particle
Click while space key pressed: create spring-connected particles
Esc key: exit

7.1. Try the Complete Version As a User¶

The dynamic simulations in this chapter deal with particles, springs, collisions between particles, and collisions between particles and walls. We will make them available as spring_mass module.

We’ll talk about how to implement the spring_mass module later, but first we’ll start by trying it out as a user.

Create a folder named spring_mass in your projects folder, and open it with VSCode. Create a file named spring_mass.py in the folder, copy and paste the code from spring_mass.py (ver 18.0) at the end of this chapter, and save it. At this point, the title bar of VSCode should show spring_mass.py - spring_mass - VSCodium (Computer Seminar I).

Then, create a simple_main.py file which will use the spring_mass module and enter the following code.

simple_main.py (ver 1.0)¶

import pygame
import spring_mass as spm


def main():
    pygame.init()
    width, height = 600, 400
    screen = pygame.display.set_mode((width, height))
    clock = pygame.time.Clock()

    world = spm.World((width, height), dt=1.0, gravity_acc=(0, 0.5))
    actor_list = []


    while True:
        frames_per_second = 60
        clock.tick(frames_per_second)

        should_quit = False
        for event in pygame.event.get():
            if event.type == pygame.QUIT:
                should_quit = True
            elif event.type == pygame.KEYDOWN and event.key == pygame.K_ESCAPE:
                should_quit = True
        if should_quit:
            break

        for a in actor_list:
            a.update()
        actor_list[:] = [a for a in actor_list if a.is_alive]

        screen.fill(pygame.Color("black"))
        for a in actor_list:
            a.draw(screen)
        pygame.display.update()

    pygame.quit()


if __name__ == "__main__":
    main()

The basic structure is the same as in particle.py ver 1.0.

The second line imports the spring_mass module. It has a long name, so I’ve made it available as spm.

The rest of the program is just the main function, an empty framework as usual with an object of type World and an actor_list (lines 11-12).

The role of World is the same as before, but the initialization arguments are slightly different. The first argument is a tuple of the width and height of the screen. The second argument is the unit time dt, and the third argument is a tuple of gravitational acceleration vectors (previously, only the y component was passed).

In this chapter, we refer to the components such as particles and springs collectively as actors, and the name “actor_list” is a generalization of “particle_list” in the previous chapters.

In lines 28-30, the actors in the actor_list are updated, and the ones that are not is_alive are deleted. In lines 33-34, each actor is drawn.

If you see a black screen and are able to press Esc to exit, you’re good to go.

Initialize the Git repository in the Source Control sidebar, and commit simple_main.py. You may or may not want to commit spring_mass.py as well.

7.1.1. Creating a Particle¶

Let’s create a single particle. In this chapter, in addition to the attributes such as position, velocity, and radius, we add mechanical properties such as mass, viscous resistance coefficient, and restitution coefficient to realize dynamic interaction. These are set by passing them as arguments during initialization. The class name has been changed from Particle to PointMass to emphasize that it has mechanical properties, which is a big change from the previous chapters.

A PointMass type object also has a reference to the World type object, just as Particle did. One major change is that the drawing process can be replaced through a Strategy pattern. Our spring_mass module defines a class called CircleDrawer, which is a class of objects that draw a circle. It is initialized by specifying the color and line width (0 for fill).

By making the drawing process replacable, it is possible to draw not only circles, but also more elaborate graphics or to paste images. Examples will be shown later.

simple_main.py (ver 2.0)¶

    world = spm.World((width, height), dt=1.0, gravity_acc=(0, 0.5))
    actor_list = []

    pos, vel = (100, 200), (10, -15)
    radius, mass = 10, 10
    viscous_damping = 0.01
    restitution = 0.95
    circle_drawer = spm.CircleDrawer("green", 0)
    actor_list.append(spm.PointMass(pos, vel, world, radius, mass,
                                    viscous_damping, restitution, circle_drawer))

    while True:

In line 18, we create a filling-with-green circle_drawer and pass it to the PointMass object initializer in line 19. The created PointMass object is appended to the actor_list.

Note that the arguments after radius have default values. See the class definition in spring_mass.py for the specific values.

7.1.2. Creating Boundaries¶

In particle.py, bouncing was the responsibility of the Particle itself (or more precisely, of the strategy object that the Particle holds as an attribute), but in this chapter, we will create a wall or a floor (a Boundary type object) that exerts a force on the particle, and the particle will receive the force and consequently bounce back.

In creating a Boundary object, the first argument to the initializer is a normal vector, the second is the coordinates of a point included in the boundary line, the third is a World object, and the fourth is a list containing the particles to be collided (any non-particles will be ignored). The normal vector is the direction to go into the wall or floor.

simple_main.py (ver 3.0)¶

    actor_list.append(spm.Boundary((1, 0), (width, 0), world, actor_list))
    actor_list.append(spm.Boundary((0, 1), (0, height), world, actor_list))

    while True:

Line 22 creates a boundary with a rightward normal through a point (width, 0), which is the wall on the right side of the screen. Line 23 creates a floor. The left wall and ceiling can be created in the same way.

If you put them in the actor_list, particles will bounce off the corresponding boundaries; if you remove them from the actor_list, particles will not bounce off and will fly away. Please try it out.

Can we create a wall that is not at the edge of the window?

This problem can be solved by, for example, avoiding exerting force on particles that have penetrated more than a certain distance threshold into the wall. However, if the threshold is not adjusted appropriately, you may experience “slipping through” depending on the speed of the particle.

Also note that the wall we create here has a “front” and a “back”. If you want a wall that can collide from both the front and the back, you need to place two walls, one facing the front and one facing the back (each with the opposite normal direction), at the same position.

7.1.3. Creating Springs¶

The spring is represented by the class Spring. At the initialization, two particles at both ends, a World type object, a spring constant, a natural length, and a drawing strategy object are passed.

LineDrawer is provided as a class of objects for drawing straight lines.

simple_main.py (ver 4.0)¶

    p1 = spm.PointMass((300, 100), (0, 0), world, drawer=circle_drawer)
    p2 = spm.PointMass((400, 120), (0, 0), world, drawer=circle_drawer)
    spring_const = 0.01
    natural_len = 20
    line_drawer = spm.LineDrawer("white", 3)
    sp12 = spm.Spring(p1, p2, world, spring_const, natural_len, line_drawer)
    actor_list += [p1, p2, sp12]

    while True:

In this example, we create two particles and connect them with a spring. It may be difficult to understand the behavior of the spring itself when gravity is present, so you may want to try setting the acceleration of gravity in line 11 to (0, 0).

Note that collisions between particles do not occur yet.

The FixedPointMass class, which inherits from PointMass, is also available. A minimal spring mass system of which one end is fixed can be placed as follows.

simple_main.py (ver 5.0)¶

    p3 = spm.FixedPointMass((350, 50), (0, 0), world, drawer=circle_drawer)
    p4 = spm.PointMass((450, 80), (0, 0), world, drawer=circle_drawer)
    sp34 = spm.Spring(p3, p4, world, spring_const, natural_len, line_drawer)
    actor_list += [p3, p4, sp34]

    while True:

7.1.4. Resolving Collisions between Particles¶

To enable collisions between particles, put an object of type CollisionResolver into the actor_list.

When initializing, pass a World type object and a list containing the collision targets. The collision will occur between all the particles in the list.

simple_main.py (ver 6.0)¶

    actor_list.append(spm.CollisionResolver(world, actor_list))

    while True:

7.1.5. Custom Drawing¶

Here is an example of how to create a drawing class by yourself. First of all, let’s look at the definitions of CircleDrawer and LineDrawer provided by the spring_mass module:

class CircleDrawer:
    def __init__(self, color, width):
        self.color = pygame.Color(color)
        self.width = width

    def __call__(self, screen, center, radius):
        pygame.draw.circle(screen, self.color, center, radius, self.width)


class LineDrawer:
    def __init__(self, color, width):
        self.color = pygame.Color(color)
        self.width = width

    def __call__(self, screen, pos1, pos2):
        pygame.draw.line(screen, self.color, pos1, pos2, self.width)

These mean that a drawing object passed to PointMass initializer is a callable object that receives the arguments screen, center, and radius, while one passed to Spring initialier is a callable object that receives screen and the positions pos1 and pos2 at both ends.

Following these examples, let’s define ImageDrawer, which works in place of CircleDrawer, to display an arbitrary image.

simple_main.py (ver 7.0)¶

import pygame
import spring_mass as spm


class ImageDrawer:
    def __init__(self, image):
        self.image = image

    def __call__(self, screen, center, radius):
        screen.blit(self.image, (center[0] - self.image.get_width() / 2,
                                 center[1] - self.image.get_height() / 2))

    image = pygame.image.load("../../assets/player/p1_walk03.png").convert_alpha()
    image_drawer = ImageDrawer(image)
    p5 = spm.PointMass((200, 300), (5, -5), world, radius=image.get_height() / 2,
                       mass=50, restitution=0.6, drawer=image_drawer)
    actor_list.append(p5)

    while True:

The definition of ImageDrawer is as shown in lines 5-11. It receives an image (a Surface object) at initialization and pastes it with the blit method each time it is __call__’ed.

In lines 49-50, it reads the player image we used in hello_pygame.py and creates an ImageDrawer type object. The rest of the code is the same as CircleDrawer.

Note that no matter what the shape of the image is, the collision calculation is done as a circle.

I would like to calculate collisions for a variaty of shapes.

This textbook only deals with collisions between circles and between circles and lines. This is because these collisions are easy to determine. The general theory is so difficult when considering collisions of more complex shapes that a whole book could be written on the subject.

Other than circles and lines, shapes that are relatively easy to handle include axis-aligned rectangles (aa-rectangle), whose edges are aligned with the coordinate axes, and capsules, which are rectangles with a semicircle attached to each end. In cases where strictness is not required, such as in games, complex shapes can often be approximated using these shapes. If you are interested, please check them out.

pygame provides a module called sprite that performs collision detection for axis-parallel rectangles.

7.2. Slinky Drop Simulation¶

Let’s try making something that looks a little more like an application. Have you ever heard of a science experiment called Slinky Drop? You can see it in the video below.

This video ends just before the results of the experiment are revealed. If you are curious about the results, you can search for the rest of the movie. Alternatively, why not trying to simulate it yourself?

Copy the first version of simple_main.py (ver 1.0) to create slinky_drop_main.py, and add springs and masses to simulate a slinky.

slinky_drop_main.py (ver 8.0)¶

import pygame
import spring_mass as spm


def create_point_mass(pos, world, fixed=False):
    vel = (0, 0)
    radius = 10
    mass = 0.1
    viscous = 0.01
    restitution = 0.95
    if fixed:
        PointMassClass = spm.FixedPointMass
    else:
        PointMassClass = spm.PointMass
    return PointMassClass(pos, vel, world, radius, mass, viscous, restitution,
                          spm.CircleDrawer("green", width=0))


def create_spring(p1, p2, world):
    spring_const = 0.02
    natural_len = 5
    return spm.Spring(p1, p2, world, spring_const, natural_len,
                      spm.LineDrawer("white", width=2))


def main():
    pygame.init()
    width, height = 600, 400
    screen = pygame.display.set_mode((width, height))
    clock = pygame.time.Clock()

    world = spm.World((width, height), dt=1.0, gravity_acc=(0, 0.5))
    actor_list = []

    p = []
    sp = []
    spacing = 30
    p.append(create_point_mass((width/2, 0), world, fixed=True))
    for k in range(1, 15):
        p.append(create_point_mass((width/2, spacing * k), world))
        sp.append(create_spring(p[k - 1], p[k], world))
    actor_list += p + sp

    while True:
        frames_per_second = 60
        clock.tick(frames_per_second)

        should_quit = False
        for event in pygame.event.get():
            if event.type == pygame.QUIT:
                should_quit = True
            elif event.type == pygame.KEYDOWN:
                if event.key == pygame.K_ESCAPE:
                    should_quit = True
                elif event.key == pygame.K_SPACE:
                    sp[0].is_alive = False
        if should_quit:
            break

        for a in actor_list:
            a.update()
        actor_list[:] = [a for a in actor_list if a.is_alive]

        screen.fill(pygame.Color("black"))
        for a in actor_list:
            a.draw(screen)
        pygame.display.update()

    pygame.quit()


if __name__ == "__main__":
    main()

Lines 5-16, 19-23: We prepare functions to create particles and springs. By calling these functions, the main program can create objects without having to worry about setting detailed arguments. Functions that create objects in this way are called factory functions.
Lines 35-42: A fixed particle is placed at the top of the center of the screen, and connected with 15 particles in series with springs.
Lines 52-56: When the Space key is pressed, the is_alive attribute of the topmost spring is set to False, simulating a release of the hand.

So, did it work as expected?

Doesn’t sp[0].is_alive = False in line 56 violate our policy of “avoid rewriting attributes from outside the class”?

If you are serious about following the policy, you should have a method in your Spring class that sets its own is_alive attribute to False (e.g. named disappear) and call it. Or maybe every Actor should have this method.

Many of you may think “Why bother with a method that just sets is_alive to False?” In fact, many Python programmers don’t care about this and rewrite attributes directly from outside the class, which some people even think is more Pythonic.

The reason we’ve tried to stick to the “don’t rewrite attributes from outside” policy in this textbook is to get you used to the idea of letting the object take care of its own responsibilities (letting the object change its own state). Once you have a good understanding of this concept, you (or your team) can decide the rules for writing the actual code.

7.3. Interactive Projectile Simulation¶

For a more complex example, we will create an interactive simulator in which we can put particles and springs with mouse and key operations, as shown in the animation at the beginning of this chapter.

We start with an empty framework. Save it as interactive_main.py.

interactive_main.py (ver 9.0)¶

import pygame
import spring_mass as spm


class AppMain:
    def __init__(self):
        pygame.init()
        width, height = 600, 400
        self.screen = pygame.display.set_mode((width, height))

        self.world = spm.World((width, height), dt=1.0, gravity_acc=(0, 0.1))
        self.actor_list = []

    def update(self):
        for a in self.actor_list:
            a.update()
        self.actor_list[:] = [a for a in self.actor_list if a.is_alive]

    def draw(self):
        self.screen.fill(pygame.Color("black"))
        for a in self.actor_list:
            a.draw(self.screen)
        pygame.display.update()

    def run(self):
        clock = pygame.time.Clock()

        while True:
            frames_per_second = 60
            clock.tick(frames_per_second)

            should_quit = False
            for event in pygame.event.get():
                if event.type == pygame.QUIT:
                    should_quit = True
                elif event.type == pygame.KEYDOWN and event.key == pygame.K_ESCAPE:
                    should_quit = True
            if should_quit:
                break

            self.update()
            self.draw()

        pygame.quit()


if __name__ == "__main__":
    AppMain().run()

7.3.1. Particles, Walls, and Collisions¶

First, let’s implement creation of particles, boundaries, and the collision resolver. These are relatively easy to implement, as you only need to place the objects.

interactive_main.py (ver 10.0)¶

import random
import pygame
import spring_mass as spm


class ActorFactory:
    def __init__(self, world, actor_list):
        self.world = world
        self.actor_list = actor_list

    def create_point_mass(self, pos, fixed=False):
        vel = (random.uniform(-10, 10), random.uniform(-10, 0))
        mass = 10
        radius = 10
        viscous = 0.01
        restitution = 0.95
        if fixed:
            PointMassClass = spm.FixedPointMass
            color = "gray"
        else:
            PointMassClass = spm.PointMass
            color = "green"
        return PointMassClass(pos, vel, self.world, radius, mass, viscous,
                              restitution, spm.CircleDrawer(color, width=0))

    def create_collision_resolver(self):
        return spm.CollisionResolver(self.world, self.actor_list)

    def create_boundary(self, name):
        width, height = self.world.size
        geometry = {"top": ((0, -1), (0, 0)),
                    "bottom": ((0, 1), (0, height)),
                    "left": ((-1, 0), (0, 0)),
                    "right": ((1, 0), (width, 0))}
        normal, point_included = geometry[name]
        return spm.Boundary(normal, point_included, self.world, self.actor_list)


class AppMain:
    def __init__(self):
        pygame.init()
        width, height = 600, 400
        self.screen = pygame.display.set_mode((width, height))

        self.world = spm.World((width, height), dt=1.0, gravity_acc=(0, 0.1))
        self.actor_list = []
        self.factory = ActorFactory(self.world, self.actor_list)

        self.actor_list.append(self.factory.create_collision_resolver())
        self.actor_list.append(self.factory.create_boundary("top"))
        self.actor_list.append(self.factory.create_boundary("bottom"))
        self.actor_list.append(self.factory.create_boundary("left"))
        self.actor_list.append(self.factory.create_boundary("right"))

    def update(self):

            should_quit = False
            for event in pygame.event.get():
                if event.type == pygame.QUIT:
                    should_quit = True
                elif event.type == pygame.KEYDOWN and event.key == pygame.K_ESCAPE:
                    should_quit = True
                elif event.type == pygame.MOUSEBUTTONDOWN:
                    self.actor_list.append(self.factory.create_point_mass(event.pos))
            if should_quit:
                break

Lines 6-36

For convenience, let’s define a class ActorFactory that has factory functions as methods. This object will hold the world and actor_list, which can be passed internally when creating an actor. This way, the main program does not need to pass these as arguments.

In lines 29-36, the Boundary factory is defined, and since specifying the top, bottom, left, and right boundaries in terms of normal vectors and point coordinates is tedious, this method allows them to be specified by name. We create a dictionary named geometry that returns a tuple of a normal vector and point coordinates associated to the key string, namely, either of “top”, “bottom”, “left”, and “right”, which we use when creating a Boundary type object.

Lines 47-53

We create a factory object, use it to create a collision resolver object and boundary objects, and add them to the actor_list. I think the creation of boundaries has become a bit cleaner now.

Lines 79-80

When a mouse button is pressed, a particle is created and added to the actor_list. At this stage, an unfixed particle is always created regardless of which button is pressed.

In lines 18 and 21, we assign classes to the variable PointMassClass. Is this allowed?

As mentioned a few times in past Q&As, in Python, classes are treated as objects, and can be assigned to variables or passed as arguments to functions, as in line 23.

In this case, the arguments used to initialize PointMass and FixedPointMass are exactly the same, so it is simpler to write them like this than to create each of them in if and else clauses.

If we implement the factory functions as the methods of AppMain, we don’t need to pass world and actor_list. It will be more convenient, won’t it?

As we have discussed many times, the key to combating large programs is to reduce the number of things that need to be considered simultaneously. Keeping classes as small as possible contributes to reducing it.

The methods of ActorFactory implement procedures that do not change the attributes of AppMain, and that’s why they don’t need to be, and are favored not to be, within AppMain.

By contrast, procedures that update AppMain’s attributes should not be driven out of the class (remember that in this textbook, we have adopted the policy of not changing attributes from outside the class). Don’t think “Oh, I can kill two birds with one stone by appending actors to actor_list in the create_point_mass method!” just because ActorFactory holds a reference to AppMain’s actor_list attribute.

A well-known rule of thumb is to separate the classes that create objects and those that use the objects. The use of factories is an embodyment of this rule.

I don’t really see the point of having multiple factory functions in a class.

The slight advantage offered is that you don’t have to pass the world or actor_list, which you need in common, to the methods every time you call them.

This is even more useful when you want to be able to replace a set of objects you want to create with another set. In the next section, we will add a method to ActorFactory to create a spring. The parameters for springs and particles are fixed to one set in this chapter, but you may want to have different sets of spring masses (e.g., a heavy and stiff set, or a soft and long natural-length set) for different purposes. You may also want to load parameters from an external file.

In such a case, if you define each spring mass set as a derived class of ActorFactory, you can easily switch between them by changing the factory created in line 47.

7.3.2. Springs¶

All that’s left is to add a way to create springs, so that the created particles are connected in order by springs while Space key is pressed.

interactive_main.py (ver 11.0)¶

    def create_spring(self, p1, p2):
        spring_const = 0.01
        natural_len = 20
        return spm.Spring(p1, p2, self.world, spring_const, natural_len,
                          spm.LineDrawer("white", width=2))

    def create_collision_resolver(self):

class AppMain:
    def __init__(self):
        pygame.init()
        width, height = 600, 400
        self.screen = pygame.display.set_mode((width, height))

        self.world = spm.World((width, height), dt=1.0, gravity_acc=(0, 0.1))
        self.actor_list = []
        self.factory = ActorFactory(self.world, self.actor_list)

        self.actor_list.append(self.factory.create_collision_resolver())
        self.actor_list.append(self.factory.create_boundary("top"))
        self.actor_list.append(self.factory.create_boundary("bottom"))
        self.actor_list.append(self.factory.create_boundary("left"))
        self.actor_list.append(self.factory.create_boundary("right"))

        self.point_mass_prev = None

    def add_connected_point_mass(self, pos, button):
        if button == 1:
            fixed = False
        elif button == 3:
            fixed = True
        else:
            return
        p = self.factory.create_point_mass(pos, fixed)
        self.actor_list.append(p)

        if self.point_mass_prev is not None:
            sp = self.factory.create_spring(p, self.point_mass_prev)
            self.actor_list.append(sp)
        if pygame.key.get_pressed()[pygame.K_SPACE]:
            self.point_mass_prev = p

    def update(self):

            should_quit = False
            for event in pygame.event.get():
                if event.type == pygame.QUIT:
                    should_quit = True
                elif event.type == pygame.KEYDOWN and event.key == pygame.K_ESCAPE:
                    should_quit = True
                elif event.type == pygame.KEYUP and event.key == pygame.K_SPACE:
                    self.point_mass_prev = None
                elif event.type == pygame.MOUSEBUTTONDOWN:
                    self.add_connected_point_mass(event.pos, event.button)
            if should_quit:
                break

Lines 26-30

We add a method to the factory to create Spring.

Line 61

In order to connect the particles in order, we need to remember the previous particle we created. It may seem a bit naive, but we’ll remember it as an AppMain attribute point_mass_prev. The default value is None (i.e., nothing is remembered).

Lines 63-77

The process when the mouse button is pressed is a bit complicated, so we’ll implement it as another method.

First, in lines 64-71, depending on the button number, unfixed or fixed particles are created and placed in actor_list.

Next, in lines 73-75, if we remember the previous particle, we connect it to the one we just created with a spring.

Finally, if the Space key is pressed, it remembers the particle created just now as the point_mass_prev attribute.

Lines 103-104

Don’t forget to reset point_mass_prev to None when you release the Space key. Can you imagine what will happen if you don’t do this? Comment out these two lines and check your answer.

Line 106

When the mouse button is pressed, the event is handled by the method add_connected_point_mass.

It is a little hard to understand the code because processing of Space key is split between run and add_connected_point_mass.

Not only at this point, but throughout this textbook, key/mouse input processing has not been abstracted enough. This is because I thought it would be difficult to learn if the structure is too elaborate.

A better design could be to keep the key/mouse input processing independent from AppMain, and let it call AppMain methods or accept state queries from AppMain as needed.

A nested series of if-elif-elif … is also not very readable code. There are many ways to improve the code, such as preparing event handling classes for each type of event and putting them in a dictionary with event.type as a key, or using a design pattern called “Chain of Responsibility”.

7.4. The spring_mass Internals¶

This is the end of the application-side examples, and we will now look into the spring_mass module, which has been given to us without explanations.

7.4.1. Basic Design¶

The basic structure of an actor is that it has methods update and draw, and an is_alive attribute.

When update is called, the actor calls the receive_force method of target objects to which it should exert force. It passes the force vector as an argument. A PointMass object is assumed to be the target of the force, and it allows reading of the attributes pos, vel, radius, mass, viscous_damping, and restitution.

The update method of PointMass calculates the force it exerts on itself, and also calculates its own motion due to the received force.

By defining each actor under this rule, the main program can update them like:

for a in actor_list:
    a.update()
actor_list[:] = [a for a in actor_list if a.is_alive]

and all the actors will work together regardless of the specific type of a (PointMass, Spring, etc.). Also:

for a in actor_list:
    a.draw(screen)

will do the necessary drawing regardless of the concrete type of a. It is the responsibility of each actor to decide whether to generate spring restoring force, restitution force, or motion calculation by update, and it is also the responsibility of each actor to decide what will be drawn by draw. Recall that this property is called polymorphism.

In Python, we can put PointMass and Spring into actor_list without any problem, but in Java and C++, where we have to specify the type for the list, I don’t think this is possible.

When we did a similar thing for Particle and ConfinedParticle in Chapter 6, we used the list of the base class Particle because ConfinedParticle was a derived class of Particle. This time, however, there is no inheritance relationship between PointMass and Spring. What should we do in this case?

What we call Actor in this chapter is actually a conceptual representation of such a common base class.

We do not create an instance of Actor itself. It is prepared for the sole purpose of defining derived classes. Such a class is called an abstract base class (ABC).

As the questioner mentioned, there is no need to define an abstract base class in Python. However, defining one makes the meaning of the program clearer and allows for type checking. If you are interested, please check out the abc module.

https://docs.python.org/ja/3/library/abc.html

Note that here:

actor_list[:] = [a for a in actor_list if a.is_alive]

[:] within the left hand side expression is mandatory. If you write actor_list = ..., actor_list would point to a different list than the one passed to Boundary and CollisionResolver objects.

7.4.2. PointMass Class¶

Start by deleting the entire contents of spring_mass.py. In the empty file, put the following.

spring_mass.py (ver 12.0)¶

import pygame


PgVector = pygame.math.Vector2


class World:
    def __init__(self, size, dt, gravity_acc):
        self.size = size
        self.dt = dt
        self.gravity_acc = PgVector(gravity_acc)


class CircleDrawer:
    def __init__(self, color, width):
        self.color = pygame.Color(color)
        self.width = width

    def __call__(self, screen, center, radius):
        pygame.draw.circle(screen, self.color, center, radius, self.width)


class LineDrawer:
    def __init__(self, color, width):
        self.color = pygame.Color(color)
        self.width = width

    def __call__(self, screen, pos1, pos2):
        pygame.draw.line(screen, self.color, pos1, pos2, self.width)

We will intensively use pygame.math.Vector2 in this module, so it is made available as PgVector in line 4 (a bit unnatural, but it is actually a stepping stone to the next chapter).

The class World doesn’t need to be explained, just note that it is slightly different from World in particle.py.

CircleDrawer and LineDrawer have already been explained.

Then, append the following code to define the PointMass class.

spring_mass.py (ver 12.0)¶

def compute_gravity_force(mass, gravity_acc):
    return mass * gravity_acc


def compute_viscous_damping_force(viscous_damping, vel):
    return -viscous_damping * vel


def integrate_symplectic(pos, vel, force, mass, dt):
    vel_new = vel + force / mass * dt
    pos_new = pos + vel_new * dt
    return pos_new, vel_new


class PointMass:
    def __init__(self, pos, vel, world, radius=10.0, mass=10.0,
                 viscous_damping=0.01, restitution=0.95, drawer=None):
        self.is_alive = True
        self.world = world
        self.drawer = drawer or CircleDrawer("blue", 1)

        self.pos = PgVector(pos)
        self.vel = PgVector(vel)
        self.radius = radius
        self.mass = mass
        self.viscous_damping = viscous_damping
        self.restitution = restitution

        self.total_force = PgVector((0, 0))

    def update(self):
        self.generate_force()
        self.move()
        self.total_force = PgVector((0, 0))

    def draw(self, screen):
        self.drawer(screen, self.pos, self.radius)

    def receive_force(self, force):
        self.total_force += PgVector(force)

    def generate_force(self):
        force_g = compute_gravity_force(self.mass, self.world.gravity_acc)
        force_v = compute_viscous_damping_force(self.viscous_damping, self.vel)
        self.receive_force(force_g + force_v)

    def move(self):
        self.pos, self.vel = \
            integrate_symplectic(self.pos, self.vel, self.total_force, self.mass, self.world.dt)

PointMass has two roles: one is to exert gravitational and viscous resistance forces on itself, and the other is to update its position and velocity according to the received forces. The former is done by calling the generate_force method from within the update method (lines 73-76).

\[\begin{split}\boldsymbol{f}_\text{gravity} &= m \boldsymbol{g}\\ \boldsymbol{f}_\text{viscous} &= - d \boldsymbol{v}\end{split}\]

where \(m\), \(d\), and \(\boldsymbol{v}\) are the mass, viscous resistance coefficient, and velocity, respectively. The calculated force vector is passed to its own receive_force method.

For the latter moving role, it has an attribute called total_force to store the total force, which is initially set to zero vector (line 60). The force vector received through the receive_force method is added to this vector (line 71). In the update method, the motion is calculated from the total force accumulated up to that time by calling the move method (lines 78-80). It is the same as particle.py except that the acceleration is calculated from the force.

\[\begin{split}\boldsymbol{a} &= \frac{\boldsymbol{f}_\text{total}}{m}\\ \Delta \boldsymbol{v} &= \boldsymbol{a} \Delta t\\ \Delta \boldsymbol{x} &= \boldsymbol{v} \Delta t\end{split}\]

Once the motion is determined, reset the attribute total_force to zero (line 65).

Separating the force calculation from the motion calculation via receive_force may seem tedious. However, by doing this, we can handle the force received from other actors such as springs and walls in a unified manner.

The remaining method is draw. It calls the drawer strategy self.drawer prepared at initialization.

Many of you may be unfamiliar with an expression like drawer or CircleDrawer("blue", 1) in line 51 where self.drawer is initialized. This means that if drawer is not None, it takes that value, and if None, it takes the value next to “or”. This is an idiom commonly used to define a value only when it is undefined.

To see how it works so far, comment out the part of interactive_main.py that creates the CollisionResolver and Boundary objects.

interactive_main.py (ver 12.0)¶

class AppMain:
    def __init__(self):
        pygame.init()
        width, height = 600, 400
        self.screen = pygame.display.set_mode((width, height))

        self.world = spm.World((width, height), dt=1.0, gravity_acc=(0, 0.1))
        self.actor_list = []
        self.factory = ActorFactory(self.world, self.actor_list)

        # self.actor_list.append(self.factory.create_collision_resolver())
        # self.actor_list.append(self.factory.create_boundary("top"))
        # self.actor_list.append(self.factory.create_boundary("bottom"))
        # self.actor_list.append(self.factory.create_boundary("left"))
        # self.actor_list.append(self.factory.create_boundary("right"))

        self.point_mass_prev = None

You should now be able to just left-click to create a particle. You may want to adjust the gravitational acceleration and initial velocity as you like.

Note that right-clicking (to create a fixed particle) or holding down space (to create a spring) will result in an error. This is because they have not been implemented yet.

7.4.3. FixedPointMass Class¶

Next, we define FixedPointMass by inheriting from PointMass. Append the following code to spring_mass.py.

spring_mass.py (ver 13.0)¶

class FixedPointMass(PointMass):
    def __init__(self, pos, vel, world, radius=10, mass=10,
                 viscous_damping=0.01, restitution=0.95, drawer=None):
        super().__init__(pos, vel, world, radius, mass,
                         viscous_damping, restitution, drawer)
        self.vel, self.mass = PgVector((0, 0)), 1e9

    def move(self):
        pass

It may seem that we only need to overwrite the move method (lines 90-91) to prevent it from moving, but actually we need to override the velocity and the mass as well, otherwise the collision calculation will not work correctly.

Copying and pasting the code of PointMass’s __init__ method is wasteful. Instead, we call PointMass’s __init__ method and reassign a zero vector to self.vel and a very large value to self.mass.

super(). __init__ is Python’s notation that refers to the __init__ method of its super class.

7.4.4. Spring Class¶

Then, let’s define Spring.

spring_mass.py (ver 14.0)¶

def compute_restoring_force(pos1, pos2, spring_const, natural_len):
    if pos1 == pos2:
        return None
    vector12 = pos2 - pos1
    distance = vector12.magnitude()
    unit_vector12 = vector12 / distance
    f1 = unit_vector12 * spring_const * (distance - natural_len)
    return f1


class Spring:
    def __init__(self, point_mass1, point_mass2, world,
                 spring_const=0.01, natural_len=0.0, drawer=None):
        self.is_alive = True
        self.world = world
        self.drawer = drawer or LineDrawer("blue", 1)

        self.p1 = point_mass1
        self.p2 = point_mass2
        self.spring_const = spring_const
        self.natural_len = natural_len

    def update(self):
        if not (self.p1.is_alive and self.p2.is_alive):
            self.is_alive = False
            return
        self.generate_force()

    def draw(self, screen):
        self.drawer(screen, self.p1.pos, self.p2.pos)

    def generate_force(self):
        f1 = compute_restoring_force(self.p1.pos, self.p2.pos, self.spring_const, self.natural_len)
        if f1 is None:
            return
        self.p1.receive_force(f1)
        self.p2.receive_force(-f1)

In the update method of the Spring class, it first checks for the existence of particles at both ends and sets its own is_alive to False if one of them does not exist (lines 117-119). If both are present, the generate_force method is called to apply the force.

The calculation of the spring restoring force itself is not so difficult. The force received by particle 1 is:

\[\boldsymbol{f}_1 = k (L - L_\text{natural}) \boldsymbol{n}_\text{1,2}\]

where \(\boldsymbol{n}_\text{1,2}\) is the unit vector from particle 1 to 2, and \(k\), \(L\), \(L_\text{natural}\) are the spring constant, the current length, and the natural length, respectively.

Lines 95-100 are basically the same as in this equation, but note that vector12 is divided by the distance to get the unit vector (line 99), and if the distance is zero, a division by zero will occur. This situation should not occur if the particles have finite size (putting aside that a point mass should not have size) and collision detection is enabled, but it can occur in the simulation. We will simply assume that no force is generated in this case.

Run interactive_main.py and hold down the Space key while placing a particle to see how the spring behaves. It may be easier to understand if you set the acceleration of gravity and initial velocity to zero.

How should we write unit tests for this?

Create a file test_spring_mass.py, write the following, and run Run All Tests from the Explorer sidebar. In this folder, pytest is not yet enabled, so you will get a popup asking you to enable it, just like in the particle folder.

Here we define a function to test generate_restoring_force. Since it is easy to prepare a “correct result” for the spring restoring force, we give it as the argument f1_expected. If None is expected (i.e., the restoring force is not calculated), the argument expects_none should be set to True. These arguments are given in @pytest.mark.parametrize.

Only the springs extending in the x and y directions are tested here. We should have tested the diagonal direction as well, but since the function under test is calculated in vector notation, it would be okay if the two independent directions moved correctly, so I omitted it. If you are concerned about this, you may want to add tests for easy-to-understand directions such as 45 degrees and 60 degrees.

test_spring_mass.py (ver 14.0)¶

import pytest
import spring_mass as spm
from spring_mass import PgVector


@pytest.mark.parametrize("pos1v, pos2v, k, nl, f1_expected, expects_none", [
    ((0, 0), (10, 0), 2, 0,  (2 * 10, 0), False),
    ((10, 0), (0, 0), 2, 0,  (-2 * 10, 0), False),
    ((0, 0), (10, 0), 2, 7,  (2 * (10 - 7), 0), False),
    ((10, 0), (0, 0), 2, 7,  (-2 * (10 - 7), 0), False),
    ((0, 0), (10, 0), 2, 13,  (-2 * (13 - 10), 0), False),
    ((10, 0), (0, 0), 2, 13,  (2 * (13 - 10), 0), False),

    ((0, 0), (0, 10), 3, 0,  (0, 3 * 10), False),
    ((2, 10), (2, 5), 4, 0,  (0, -5 * 4), False),
    ((0, 0), (0, 10), 3, 5,  (0, 3 * 5), False),
    ((0, 0), (0, 10), 3, 12,  (0, -3 * 2), False),

    ((0, 0), (0, 0), 0, 0, (0, 0), True), # equal positions
    ((12, 5), (12, 5), 0, 0, (0, 0), True), # equal positions
])
def test_restoring_force(pos1v, pos2v, k, nl, f1_expected, expects_none):
    pos1 = PgVector(pos1v)
    pos2 = PgVector(pos2v)
    f1 = spm.compute_restoring_force(pos1, pos2, k, nl)
    if expects_none:
        assert f1 is None
    else:
        assert f1 == PgVector(f1_expected)

When I’m given a mathematical expression like this, I can understand that it will be a program like this, but if I’m asked to think of a mathematical expression by myself, it’s difficult. In particular, I can’t confidently figure out how to deal with the sign of the force depending on whether the spring length is longer or shorter than the natural length.

It may sound surprising, but there are cases where it is effective to write tests first, and then write the target program so that it passes the tests. This is called test driven development (TDD).

Let’s consider the example of the spring restoration force. First, read the explanation of unit testing in the previous Q&A if you haven’t.

Now, if the natural length is 0, you can write the code without thinking about anything difficult. So, first, we will prepare a test only for that case. Let’s consider only the case where it extends in the x direction. The expected answer f1_expected can be calculated by hand:

@pytest.mark.parametrize("pos1v, pos2v, k, nl, f1_expected, expects_none", [
    ((0, 0), (10, 0), 2, 0,  (2 * 10, 0), False),
    ((10, 0), (0, 0), 2, 0,  (-2 * 10, 0), False),
])
def test_restoring_force(pos1v, pos2v, k, nl, f1_expected, expects_none):
    pos1 = PgVector(pos1v)
    pos2 = PgVector(pos2v)
    f1 = spm.compute_restoring_force(pos1, pos2, k, nl)
    if expects_none:
        assert f1 is None
    else:
        assert f1 == PgVector(f1_expected)

And it is easy to write compute_restoring_force that will pass this test:

def compute_restoring_force(pos1, pos2, spring_const, natural_len):
    vector12 = pos2 - pos1
    distance = vector12.magnitude()
    unit_vector12 = vector12 / distance
    f1 = unit_vector12 * spring_const * distance
    return f1

Next, let’s add a test for the case where the natural length is non-zero:

@pytest.mark.parametrize("pos1v, pos2v, k, nl, f1_expected, expects_none", [
    ((0, 0), (10, 0), 2, 0,  (2 * 10, 0), False),
    ((10, 0), (0, 0), 2, 0,  (-2 * 10, 0), False),
    ((0, 0), (10, 0), 2, 7,  (2 * (10 - 7), 0), False),
    ((10, 0), (0, 0), 2, 7,  (-2 * (10 - 7), 0), False),

We have prepared the case where the natural length nl is shorter than the spring length = 10. The sign of the force generated in these cases is easily known, so we specify it as f1_expected.

The compute_restoring_force we previously defined should not pass the test for the new two cases. It is not difficult to fix it to pass:

def compute_restoring_force(pos1, pos2, spring_const, natural_len):
    vector12 = pos2 - pos1
    distance = vector12.magnitude()
    unit_vector12 = vector12 / distance
    f1 = unit_vector12 * spring_const * (distance - natural_len)
    return f1

In a similar way, let’s add a test for the case where the natural length is longer than the spring length:

@pytest.mark.parametrize("pos1v, pos2v, k, nl, f1_expected, expects_none", [
    ((0, 0), (10, 0), 2, 0,  (2 * 10, 0), False),
    ((10, 0), (0, 0), 2, 0,  (-2 * 10, 0), False),
    ((0, 0), (10, 0), 2, 7,  (2 * (10 - 7), 0), False),
    ((10, 0), (0, 0), 2, 7,  (-2 * (10 - 7), 0), False),
    ((0, 0), (10, 0), 2, 13,  (-2 * (13 - 10), 0), False),
    ((10, 0), (0, 0), 2, 13,  (2 * (13 - 10), 0), False),

When you run the tests, these new test cases should also pass. It looks like the previous code is fine.

Let’s add a test for the y direction. Let’s also try slightly different positions:

@pytest.mark.parametrize("pos1v, pos2v, k, nl, f1_expected, expects_none", [
    ((0, 0), (10, 0), 2, 0,  (2 * 10, 0), False),
    ((10, 0), (0, 0), 2, 0,  (-2 * 10, 0), False),
    ((0, 0), (10, 0), 2, 7,  (2 * (10 - 7), 0), False),
    ((10, 0), (0, 0), 2, 7,  (-2 * (10 - 7), 0), False),
    ((0, 0), (10, 0), 2, 13,  (-2 * (13 - 10), 0), False),
    ((10, 0), (0, 0), 2, 13,  (2 * (13 - 10), 0), False),

    ((0, 0), (0, 10), 3, 0,  (0, 3 * 10), False),
    ((2, 10), (2, 5), 4, 0,  (0, -5 * 4), False),
    ((0, 0), (0, 10), 3, 5,  (0, 3 * 5), False),
    ((0, 0), (0, 10), 3, 12,  (0, -3 * 2), False),

It seems to pass without any problems.

At this point, we can feel confident that the basic logic seems to be okay. If you are concerned, you can add a diagonal test.

The next step is to check if there are any exceptional conditions. If you reread the compute_restoring_force you wrote, you will see that it executes division by distance, which may cause division by zero. Let’s add it as a test case:

@pytest.mark.parametrize("pos1v, pos2v, k, nl, f1_expected, expects_none", [
    ((0, 0), (10, 0), 2, 0,  (2 * 10, 0), False),
    ((10, 0), (0, 0), 2, 0,  (-2 * 10, 0), False),
    ((0, 0), (10, 0), 2, 7,  (2 * (10 - 7), 0), False),
    ((10, 0), (0, 0), 2, 7,  (-2 * (10 - 7), 0), False),
    ((0, 0), (10, 0), 2, 13,  (-2 * (13 - 10), 0), False),
    ((10, 0), (0, 0), 2, 13,  (2 * (13 - 10), 0), False),

    ((0, 0), (0, 10), 3, 0,  (0, 3 * 10), False),
    ((2, 10), (2, 5), 4, 0,  (0, -5 * 4), False),
    ((0, 0), (0, 10), 3, 5,  (0, 3 * 5), False),
    ((0, 0), (0, 10), 3, 12,  (0, -3 * 2), False),

    ((0, 0), (0, 0), 0, 0, (0, 0), True), # equal positions
    ((12, 5), (12, 5), 0, 0, (0, 0), True), # equal positions

These test cases do not pass. Fix it to return None when the distance is 0:

def compute_restoring_force(pos1, pos2, spring_const, natural_len):
    if pos1 == pos2:
        return None
    vector12 = pos2 - pos1
    distance = vector12.magnitude()
    unit_vector12 = vector12 / distance
    f1 = unit_vector12 * spring_const * (distance - natural_len)
    return f1

Now we have reached to the ver 14.0 code.

The knack for developing in this way is to implement the core of the problem straightforwardly at first, ignoring exceptional cases or cases that require detailed consideration, and then refine it by adding tests for the finer cases.

7.4.5. FragileSpring Class¶

As an example of extending Spring, let’s define FragileSpring, a spring that breaks when a force exceeding a certain threshold is applied.

spring_mass.py (ver 15.0)¶

class FragileSpring(Spring):
    def __init__(self, point_mass1, point_mass2, world,
                 spring_const=0.01, natural_len=0.0, drawer=None,
                 break_threshold=1e9):
        super().__init__(point_mass1, point_mass2, world, spring_const,
                         natural_len, drawer)
        self.break_threshold = break_threshold

    def generate_force(self):
        f1 = compute_restoring_force(self.p1.pos, self.p2.pos, self.spring_const, self.natural_len)
        if f1 is None:
            return
        self.p1.receive_force(f1)
        self.p2.receive_force(-f1)
        if f1.magnitude() > self.break_threshold:
            self.is_alive = False

In the __init__ method of FragileSpring, we call __init__ of the super class and add the break_threshold attribute (lines 137-139).

In the generate_force method, we would like to call the super class generate_force… but we don’t have a way to get the force calculation results, so we are forced to copy-paste the same code. Then add the rupture process (lines 147-148).

In interactive_main.py, modify the factory’s create_spring method to create FragileSpring instead of Spring.

interactive_main.py (ver 15.0)¶

    def create_spring(self, p1, p2):
        spring_const = 0.01
        natural_len = 20
        break_threshold = 5.0
        return spm.FragileSpring(p1, p2, self.world, spring_const, natural_len,
                                 spm.LineDrawer("white", width=2), break_threshold)

With these constant setting, if you make a spring that is as long as the edge of the screen, it should disappear.

7.4.6. CollisionResolver Class¶

Next is CollisionResolver. It may sound difficult, but the basic idea is the same as for other actors, and the responsibility is to call the receive_force method of the actor to which the force should be applied and pass the force vector.

spring_mass.py (ver 16.0)¶

def is_point_mass(actor):
    return isinstance(actor, PointMass)


def compute_impact_force_between_points(p1, p2, dt):
    if (p1.pos - p2.pos).magnitude() > p1.radius + p2.radius:
        return None
    if p1.pos == p2.pos:
        return None
    normal = (p2.pos - p1.pos).normalize()
    v1 = p1.vel.dot(normal)
    v2 = p2.vel.dot(normal)
    if v1 < v2:
        return None
    e = p1.restitution * p2.restitution
    m1, m2 = p1.mass, p2.mass
    f1 = normal * (-(e + 1) * v1 + (e + 1) * v2) / (1/m1 + 1/m2) / dt
    return f1


class CollisionResolver:
    def __init__(self, world, actor_list, target_condition=None, drawer=None):
        self.is_alive = True
        self.world = world
        self.drawer = drawer or (lambda surface: None)

        self.actor_list = actor_list
        if target_condition is None:
            self.target_condition = is_point_mass
        else:
            self.target_condition = target_condition

    def update(self):
        self.generate_force()

    def draw(self, surface):
        self.drawer(surface)


    def generate_force(self):
        plist = [a for a in self.actor_list if self.target_condition(a)]
        n = len(plist)
        for i in range(n):
            for j in range(i + 1, n):
                p1, p2 = plist[i], plist[j]
                f1 = compute_impact_force_between_points(p1, p2, self.world.dt)
                if f1 is None:
                    continue
                p1.receive_force(f1)
                p2.receive_force(-f1)

Line 175 initializes the drawing strategy self.drawer. If the argument drawer is not None, it takes that value. If None, it takes a function that draws nothing written by a lambda expression.

The update method simply calls the generate_force method (line 184). The generate_force method first lists the actors to which the force should be applied (line 191), determines the collision between them, and generates the force if there is a collision.

In line 191, only the actors in the actor_list attribute that satisfy the target condition are extracted to the plist variable. The condition is determined by the target_condition function, which by default (lines 178-179) determines if the actor is an instance of the PointMass class (lines 151-152). The condition is given in the Strategy pattern, so it can be replaced as needed.

The isinstance(x, A) function used here is provided by Python and is true if the object x is an instance of class A. However, its feature is that it becomes also true for instances of A’s sub classes, not just A itself. Therefore, FixedPointMass is also subject to collision calculation.

After listing the n target objects in the plist, it lists up the combinations that take out two different objects from the list. Say one is the i-th object in the plist and the other is the j-th object, then you can get all combinations without duplicates by enumerating j from i+1 to n-1 while enumerating i from 0 to n-1 (lines 193-194).

The computation of collision is done by compute_impact_force_between_points function (lines 155-168). First, if the distance between the centers is greater than the sum of the radii, there is no collision and we exclude it (lines 156-157). If there is a collision, we assume that the elastic collision force is generated along the vector connecting the centers.

Now, it is time for high school physics. Consider a one-dimensional motion with only a component in the direction from particle 1 to particle 2. If the masses and velocities of particles 1 and 2 are \(m_1\), \(m_2\), \(v_1\), and \(v_2\), respectively, the momentum conservation law is

\[m_1 v_1 + m_2 v_2 = m_1 (v_1 + \Delta v_1) + m_2 (v_2 + \Delta v_2)\]

and with the restitution coefficient e,

\[(v_1 + \Delta v_1) - (v_2 + \Delta v_2) = - e (v_1 - v_2)\]

is required to be satisfied. The force that must be applied to particles 1 and 2 to induce the velocity change that satisfies these conditions is given by

\[\begin{split}f_1 &= m_1 \frac{\Delta v_1}{\Delta t}\\ f_2 &= - f_1\end{split}\]

By solving these equations, we get the following force:

\[f_1 = \frac{-(e + 1) v_1 + (e + 1) v_2}{(\frac{1}{m_1} + \frac{1}{m_2}) \Delta t}\]

So, according to this equation, we will generate a force in the direction from particle 1 to 2. Find the unit vector in this direction (line 160), and calculate the velocity component in that direction as the inner product using the dot method of pygame.math.Vector2 (lines 161-162). Then multiply this unit direction vector by the magnitude of the force and you’re good to go (line 167).

However, there are a few details that should be considered.

First of all, the direction of the force cannot be determined if the positions of particles 1 and 2 are exactly the same. As in the case of the spring, this is also an improbable situation, but it can happen in simulations. We again simply ignore it (lines 158-159).
Next, we need to deal with the phenomenon of multiple collisions, which occurred when we considered collisions with walls in particle.py. As in particle.py, we consider the velocity condition, i.e., we exclude cases where the relative velocity of the two objects is moving away from each other (lines 163-164).
Finally, we need to know how to give the restitution coefficient. Essentially, the restitution coefficient must be given for each pair of colliding bodies, and cannot be defined for a single particle. However, since it is troublesome for users to define the restitution coefficient for each pair of colliding bodies, many physics simulators employ a convenient approach of defining the restitution coefficient (or something like it) for each individual body, and then using some rule to calculate the restitution coefficient for a pair. Here we use the product of the restitution attributes of both. This is not physically grounded at all.

When you are done, reactivate the CollisionResolver in interactive_main.py and check the behavior.

interactive_main.py (ver 16.0)¶

        self.factory = ActorFactory(self.world, self.actor_list)

        self.actor_list.append(self.factory.create_collision_resolver())
        # self.actor_list.append(self.factory.create_boundary("top"))
        # self.actor_list.append(self.factory.create_boundary("bottom"))

I was wondering the same thing about the spring, but why does it return None when no force is generated? Wouldn’t the code be cleaner if it returned PgVector((0, 0))?

In the case of springs it’s not so much a problem, but in the case of collisions there can be many combinations to check (proportional to the square of the number of point masses). Returning PgVector((0, 0)) will give the same result, but it will be slower because the receive_force method will be called unnecessarily often.

What kind of test should we write?

Since the equations are more complicated than the spring restoring force, it is tedious to give each “correct” result from outside. However, if we write the equation computed in the function under test, compute_impact_force_between_points, also in the test function as it is, the significance of the test will be weak.

We will change our mind and calculate the velocities before and after the update by the force to confirm that momentum is conserved and that the elastic collision relation is valid.

test_spring_mass.py (ver 16.0)¶

@pytest.mark.parametrize("p1v, p2v, expects_none", [
    (((0, 0), (5, 0), 20, 10, 1.0),   # (pos, vel, r, m, e)
     ((20, 0), (-5, 0), 20, 10, 1.0), False),
    (((0, 0), (5, 0), 20, 15, 0.2),
     ((20, 0), (-5, 0), 20, 30, 0.8), False),

    (((0, 0), (5, 0), 20, 10, 1.0),
     ((120, 0), (-5, 0), 20, 10, 1.0), True), # not in contact
    (((0, 0), (5, 0), 20, 10, 1.0),
     ((0, 0), (-5, 0), 20, 10, 1.0), True), # equal positions
    (((0, 0), (-5, 0), 20, 10, 1.0),
     ((20, 0), (5, 0), 20, 10, 1.0), True), # moving apart
])
def test_impact_force_between_points(p1v, p2v, expects_none):
    dt, g = 1.0, 0
    w = spm.World((600, 400), dt, g)
    pos1, vel1, r1, m1, e1 = p1v
    pos2, vel2, r2, m2, e2 = p2v

    p1 = spm.PointMass(pos1, vel1, w,
                radius=r1, mass=m1, restitution=e1)
    p2 = spm.PointMass(pos2, vel2, w,
                radius=r2, mass=m2, restitution=e2)
    f1 = spm.compute_impact_force_between_points(p1, p2, dt)
    if expects_none:
        assert f1 is None
        return

    vel1_old = PgVector(vel1)
    vel2_old = PgVector(vel2)
    vel1_new = vel1_old + f1 / m1 * dt
    vel2_new = vel2_old - f1 / m2 * dt
    e = e1 * e2
    assert vel1_new - vel2_new == -e * (vel1_old - vel2_old)
    assert m1 * vel1_new + m2 * vel2_new == m1 * vel1_old + m2 * vel2_old

7.4.7. Boundary Class¶

I think you are getting tired by now, but the rest is Boundary class.

spring_mass.py (ver 17.0)¶

def compute_impact_force_by_fixture(p, normal, point_included, dt):
    invasion = normal.dot(p.pos - point_included)
    if invasion + p.radius > 0 and normal.dot(p.vel) > 0:
        e = p.restitution
        v = normal.dot(p.vel)
        m = p.mass
        f = normal * (-(e + 1) * v) * m / dt
    else:
        f = None
    return f


class Boundary:
    def __init__(self, normal, point_included, world, actor_list,
                 target_condition=None, drawer=None):
        self.is_alive = True
        self.world = world
        self.drawer = drawer or (lambda surface: None)

        self.normal = PgVector(normal).normalize()
        self.point_included = PgVector(point_included)
        self.actor_list = actor_list
        if target_condition is None:
            self.target_condition = is_point_mass
        else:
            self.target_condition = target_condition

    def update(self):
        self.generate_force()

    def draw(self, surface):
        self.drawer(surface)


    def generate_force(self):
        plist = [a for a in self.actor_list if self.target_condition(a)]
        for p in plist:
            f = compute_impact_force_by_fixture(p, self.normal, self.point_included, self.world.dt)
            if f is None:
                continue
            p.receive_force(f)

The structure is similar to that of CollisionResolver, but a little simpler. In the generate_force method, collision targets are listed (line 238), and collision is calculated for each of them.

Let the normal vector of the boundary surface denoted by \(\boldsymbol{n}\), and the coordinates of a point included in the boundary surface by \(\boldsymbol{p}_0\). Recall that the normal vector is defined to be in the direction towards the interior of the boundary.

The normal component of the invation of the center position \(\boldsymbol{x}\) into the boundary is given by the inner product \(\boldsymbol{n} \cdot (\boldsymbol{x} - \boldsymbol{p}_0)\). If this plus the offset by the radius is positive and the normal component of the velocity is positive, then a collision has occurred (lines 204-205).

We could develop the equations for the force again, but it is quicker to consider the case with \(v_2 = 0\) and the limit \(m_2 \rightarrow \infty\) in the equation we used for the collision between the particles (lines 206-209).

In interactive_main.py, you can play with building and removing walls and floors.

interactive_main.py (ver 17.0)¶

        self.actor_list.append(self.factory.create_collision_resolver())
        self.actor_list.append(self.factory.create_boundary("top"))
        self.actor_list.append(self.factory.create_boundary("bottom"))
        self.actor_list.append(self.factory.create_boundary("left"))
        self.actor_list.append(self.factory.create_boundary("right"))

What about tests?

Again, since it is not interesting to test the force returned by the function under test as it is, we will calculate the velocities before and after the update by force, and confirm that the elastic collision relationship is valid.

test_spring_mass.py (ver 17.0)¶

@pytest.mark.parametrize("pv, n, pinc, vel_expected, expects_none", [
    (((600, 200), (5, 3), 20, 10, 1.0),   # (pos, vel, r, m, e)
     (1, 0), (600, 0), (-5, 3), False),
    (((15, 200), (-5, 3), 20, 10, 0.8),
     (-1, 0), (0, 0), (0.8 * 5, 3), False),
    (((300, 400), (-5, 3), 20, 10, 0.6),
     (0, 1), (0, 400), (-5, - 0.6 * 3), False),

    (((600, 200), (-5, 3), 20, 10, 1.0),
     (1, 0), (600, 0), (0, 0), True), # not moving toward warll
    (((580, 200), (5, 3), 20, 10, 1.0),
     (1, 0), (600, 0), (0, 0), True), # not in contact
])
def test_impact_force_by_fixture(pv, n, pinc, vel_expected, expects_none):
    dt, g = 1.0, 0
    w = spm.World((600, 400), dt, g)
    pos, vel, r, m, e = pv
    p = spm.PointMass(pos, vel, w, radius=r, mass=m, restitution=e)
    normal = PgVector(n)
    point_included = PgVector(pinc)
    f = spm.compute_impact_force_by_fixture(p, normal, point_included, dt)
    if expects_none:
        assert f is None
        return

    vel_old = PgVector(vel)
    vel_new = vel_old + f / m * dt
    assert vel_new == PgVector(vel_expected)

7.4.8. Avoiding Sinking into Floor¶

There is one last problem: sinking into the floor, which was also a problem in particle.py.

As in particle.py, we want to avoid this by forcing the y-coordinate to be on the floor, but this time the problem is a bit more subtle. Is this a PointMass responsibility, or a Boundary responsibility?

It is essentially the responsibility of PointMass, since it concerns the calculation of the position of a particle. However, a PointMass object has no way of knowing when it has collided with the floor, because it receives the force from the spring, the restitution force from the other particles, and the restitution force from the floor through the receive_force method without distinguishing between them.

Then, should Boundary change the y-coordinate of a particle? I don’t want to do it because it is against the policy to avoid changing the attributes of an object from the outside.

If we want to expand the spring_mass module further, similar problems can occur in other cases. For example, let’s say we want to define FragilePointMass that breaks after a certain number of collisions with other particles. However, there is no way to distinguish the receipt of a force from a collision with another particle from the receipt of other forces.

In short, we need a way to transmit information other than force from the generate_force side to the receive_force side. Here, we will implement the following messaging mechanism.

spring_mass.py (ver 18.0)¶

import pygame


PgVector = pygame.math.Vector2


class World:
    def __init__(self, size, dt, gravity_acc):
        self.size = size
        self.dt = dt
        self.gravity_acc = PgVector(gravity_acc)


class CircleDrawer:
    def __init__(self, color, width):
        self.color = pygame.Color(color)
        self.width = width

    def __call__(self, screen, center, radius):
        pygame.draw.circle(screen, self.color, center, radius, self.width)


class LineDrawer:
    def __init__(self, color, width):
        self.color = pygame.Color(color)
        self.width = width

    def __call__(self, screen, pos1, pos2):
        pygame.draw.line(screen, self.color, pos1, pos2, self.width)


def compute_gravity_force(mass, gravity_acc):
    return mass * gravity_acc


def compute_viscous_damping_force(viscous_damping, vel):
    return -viscous_damping * vel


def integrate_symplectic(pos, vel, force, mass, dt):
    vel_new = vel + force / mass * dt
    pos_new = pos + vel_new * dt
    return pos_new, vel_new


class PointMass:
    def __init__(self, pos, vel, world, radius=10.0, mass=10.0,
                 viscous_damping=0.01, restitution=0.95, drawer=None):
        self.is_alive = True
        self.world = world
        self.drawer = drawer or CircleDrawer("blue", 1)

        self.pos = PgVector(pos)
        self.vel = PgVector(vel)
        self.radius = radius
        self.mass = mass
        self.viscous_damping = viscous_damping
        self.restitution = restitution

        self.total_force = PgVector((0, 0))
        self.message_list = []

    def update(self):
        self.generate_force()
        self.move()
        self.total_force = PgVector((0, 0))
        self.message_list.clear()

    def draw(self, screen):
        self.drawer(screen, self.pos, self.radius)

    def receive_force(self, force):
        self.total_force += PgVector(force)

    def receive_message(self, msg):
        self.message_list.append(msg)

    def generate_force(self):
        force_g = compute_gravity_force(self.mass, self.world.gravity_acc)
        force_v = compute_viscous_damping_force(self.viscous_damping, self.vel)
        self.receive_force(force_g + force_v)

    def move(self):
        self.pos, self.vel = \
            integrate_symplectic(self.pos, self.vel, self.total_force, self.mass, self.world.dt)

        for msg in self.message_list:
            if msg["type"] == "floor_hit" and self.vel.y > 0:
                # constrain y on or above floor
                self.pos.y = msg["y"] - self.radius


class FixedPointMass(PointMass):
    def __init__(self, pos, vel, world, radius=10.0, mass=10.0,
                 viscous_damping=0.01, restitution=0.95, drawer=None):
        super().__init__(pos, vel, world, radius, mass,
                         viscous_damping, restitution, drawer)
        self.vel, self.mass = PgVector((0, 0)), 1e9

    def move(self):
        pass


def compute_restoring_force(pos1, pos2, spring_const, natural_len):
    if pos1 == pos2:
        return None
    vector12 = pos2 - pos1
    distance = vector12.magnitude()
    unit_vector12 = vector12 / distance
    f1 = unit_vector12 * spring_const * (distance - natural_len)
    return f1


class Spring:
    def __init__(self, point_mass1, point_mass2, world,
                 spring_const=0.01, natural_len=0.0, drawer=None):
        self.is_alive = True
        self.world = world
        self.drawer = drawer or LineDrawer("blue", 1)

        self.p1 = point_mass1
        self.p2 = point_mass2
        self.spring_const = spring_const
        self.natural_len = natural_len

    def update(self):
        if not (self.p1.is_alive and self.p2.is_alive):
            self.is_alive = False
            return
        self.generate_force()

    def draw(self, screen):
        self.drawer(screen, self.p1.pos, self.p2.pos)

    def generate_force(self):
        f1 = compute_restoring_force(self.p1.pos, self.p2.pos, self.spring_const, self.natural_len)
        if f1 is None:
            return
        self.p1.receive_force(f1)
        self.p2.receive_force(-f1)


class FragileSpring(Spring):
    def __init__(self, point_mass1, point_mass2, world,
                 spring_const=0.01, natural_len=0.0, drawer=None,
                 break_threshold=1e9):
        super().__init__(point_mass1, point_mass2, world, spring_const,
                         natural_len, drawer)
        self.break_threshold = break_threshold

    def generate_force(self):
        f1 = compute_restoring_force(self.p1.pos, self.p2.pos, self.spring_const, self.natural_len)
        if f1 is None:
            return
        self.p1.receive_force(f1)
        self.p2.receive_force(-f1)
        if f1.magnitude() > self.break_threshold:
            self.is_alive = False


def is_point_mass(actor):
    return isinstance(actor, PointMass)


def compute_impact_force_between_points(p1, p2, dt):
    if (p1.pos - p2.pos).magnitude() > p1.radius + p2.radius:
        return None
    if p1.pos == p2.pos:
        return None
    normal = (p2.pos - p1.pos).normalize()
    v1 = p1.vel.dot(normal)
    v2 = p2.vel.dot(normal)
    if v1 < v2:
        return None
    e = p1.restitution * p2.restitution
    m1, m2 = p1.mass, p2.mass
    f1 = normal * (-(e + 1) * v1 + (e + 1) * v2) / (1/m1 + 1/m2) / dt
    return f1


class CollisionResolver:
    def __init__(self, world, actor_list, target_condition=None, drawer=None):
        self.is_alive = True
        self.world = world
        self.drawer = drawer or (lambda surface: None)

        self.actor_list = actor_list
        if target_condition is None:
            self.target_condition = is_point_mass
        else:
            self.target_condition = target_condition

    def update(self):
        self.generate_force()

    def draw(self, surface):
        self.drawer(surface)


    def generate_force(self):
        plist = [a for a in self.actor_list if self.target_condition(a)]
        n = len(plist)
        for i in range(n):
            for j in range(i + 1, n):
                p1, p2 = plist[i], plist[j]
                f1 = compute_impact_force_between_points(p1, p2, self.world.dt)
                if f1 is None:
                    continue
                p1.receive_force(f1)
                p2.receive_force(-f1)


def compute_impact_force_by_fixture(p, normal, point_included, dt):
    invasion = normal.dot(p.pos - point_included)
    if invasion + p.radius > 0 and normal.dot(p.vel) > 0:
        e = p.restitution
        v = normal.dot(p.vel)
        m = p.mass
        f = normal * (-(e + 1) * v) * m / dt
    else:
        f = None
    return f


class Boundary:
    def __init__(self, normal, point_included, world, actor_list,
                 target_condition=None, drawer=None):
        self.is_alive = True
        self.world = world
        self.drawer = drawer or (lambda surface: None)

        self.normal = PgVector(normal).normalize()
        self.point_included = PgVector(point_included)
        self.actor_list = actor_list
        if target_condition is None:
            self.target_condition = is_point_mass
        else:
            self.target_condition = target_condition

    def update(self):
        self.generate_force()

    def draw(self, surface):
        self.drawer(surface)


    def is_floor(self):
        return self.normal == PgVector((0, 1))

    def generate_force(self):
        plist = [a for a in self.actor_list if self.target_condition(a)]
        for p in plist:
            f = compute_impact_force_by_fixture(p, self.normal, self.point_included, self.world.dt)
            if f is None:
                continue
            p.receive_force(f)
            if self.is_floor():
                p.receive_message({"type": "floor_hit", "y": self.point_included.y})

Lines 61, 67: As an attribute of the force receiver, we prepare a list message_list to store messages separately from the total_force. This list is initially set to be empty (line 61) and is reset to empty after each update (line 67).
Lines 75-76: In addition to the receive_force method, we provide a receive_message method.
Lines 258-259: The method to send a message is explained first. Boundary’s generate_force method generates a force in line 256 and then sends a message to the destination if it is a floor object. The message is created as a Python dictionary with the value for the “type” key set to “floor_hit” and the value for the “y” key set to the y coordinate.
Lines 87-91: Let’s get back to PointMass. The processing of messages is similar to the that of events in pygame: it takes one message at a time from the message_list, and if its “type” is “floor_hit” and the particle has a downward velocity component, it updates the y-coordinate using the value of the “y” key in the message.

Isn’t there a more fundamental solution to this problem of sinking?

It is not easy to implement, but I will just explain the basic idea.

A particle should satisfy the constraint that it cannot enter a wall. You can define a penalty value for not satisfying this constraint (e.g., the square of the penetration), and adjust the position to minimize the penalty for each motion calculation.

Similar constraints should be satisfied not only against walls, but also between particles. (Note that the current program does not take this into account at all; you can see this by arranging several FixedPointMass’es in a mortar shape and dropping a PointMass into them)

Many of other constraints can also be expressed as distance constraints, such as when two particles are connected by a rigid link or by a non-retractable string.

When multiple constraints exist between multiple moving bodies, it is necessary to optimize the position of the entire system of particles such that the overall penalty is minimized. For example, this can be formulated as a nonlinear least-squares problem.

Compared to Spring and CollisionResolver, PointMass seems to have a quite complex structure.

PointMass is responsible for two things at once: force generation and motion calculation. From a physics point of view, there are two aspects of mass, gravitational mass and inertial mass, which are understood to be essentially different but coincide with each other.

It would be desirable to separate at least these two responsibilities into different classes (e.g., GravityGenerator and MovingPoint), and have the PointMass class keep references to each object (i.e. object composition) to divide the processing.

Another possible design is to define PointMass by inheriting from these two classes. However, inheriting multiple base classes (called multiple inheritance) generally makes the design difficult, and it is better to avoid it if possible. In some languages, such as Java, multiple inheritance is prohibited.

7.5. Exercise¶

Problem 7-1

Define your own LineDrawer such that Spring is displayed as a dotted line. For example, you can divide a line segment connecting two end points into an appropriate number of segments, and display a circle for each segment.

Problem 7-2

Consider that you have a spring with constant \(k\) connected to a fixed point and an unfixed mass \(m\) in a world without gravity. A sinusoidal force is applied to the unfixed mass to cause a forced vibration. Simulate the behavior when the angular frequency of the forced vibration is close to or far from the resonant angular frequency \(\sqrt{k / m}\).

(Hint: A simple harmonic motion with angular frequency \(\omega\) can be expressed as \(\cos(\omega n \Delta t)\), where \(n\) is the frame number and \(\Delta t\) is the time step per frame. The amplitude should be set appropriately)

Problem 7-3

Create a collision resolver class that inherits from CollisionResolver and displays the total number of collisions as a string in the upper left corner of the screen.

(Hint: For example, you can add an attribute to hold the total number of collisions, override generate_force to count it up, and override draw to display it)

Problem 7-4

By referring to the Spring class implementation, create an actor class that generates attraction or repulsion force inversely proportional to the squared distance between two particles.