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.
1import pygame
2import spring_mass as spm
3
4
5def main():
6 pygame.init()
7 width, height = 600, 400
8 screen = pygame.display.set_mode((width, height))
9 clock = pygame.time.Clock()
10
11 world = spm.World((width, height), dt=1.0, gravity_acc=(0, 0.5))
12 actor_list = []
13
14
15 while True:
16 frames_per_second = 60
17 clock.tick(frames_per_second)
18
19 should_quit = False
20 for event in pygame.event.get():
21 if event.type == pygame.QUIT:
22 should_quit = True
23 elif event.type == pygame.KEYDOWN and event.key == pygame.K_ESCAPE:
24 should_quit = True
25 if should_quit:
26 break
27
28 for a in actor_list:
29 a.update()
30 actor_list[:] = [a for a in actor_list if a.is_alive]
31
32 screen.fill(pygame.Color("black"))
33 for a in actor_list:
34 a.draw(screen)
35 pygame.display.update()
36
37 pygame.quit()
38
39
40if __name__ == "__main__":
41 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.
11 world = spm.World((width, height), dt=1.0, gravity_acc=(0, 0.5))
12 actor_list = []
13
14 pos, vel = (100, 200), (10, -15)
15 radius, mass = 10, 10
16 viscous_damping = 0.01
17 restitution = 0.95
18 circle_drawer = spm.CircleDrawer("green", 0)
19 actor_list.append(spm.PointMass(pos, vel, world, radius, mass,
20 viscous_damping, restitution, circle_drawer))
21
22 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.
21
22 actor_list.append(spm.Boundary((1, 0), (width, 0), world, actor_list))
23 actor_list.append(spm.Boundary((0, 1), (0, height), world, actor_list))
24
25 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.
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.
24
25 p1 = spm.PointMass((300, 100), (0, 0), world, drawer=circle_drawer)
26 p2 = spm.PointMass((400, 120), (0, 0), world, drawer=circle_drawer)
27 spring_const = 0.01
28 natural_len = 20
29 line_drawer = spm.LineDrawer("white", 3)
30 sp12 = spm.Spring(p1, p2, world, spring_const, natural_len, line_drawer)
31 actor_list += [p1, p2, sp12]
32
33 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.
32
33 p3 = spm.FixedPointMass((350, 50), (0, 0), world, drawer=circle_drawer)
34 p4 = spm.PointMass((450, 80), (0, 0), world, drawer=circle_drawer)
35 sp34 = spm.Spring(p3, p4, world, spring_const, natural_len, line_drawer)
36 actor_list += [p3, p4, sp34]
37
38 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.
37
38 actor_list.append(spm.CollisionResolver(world, actor_list))
39
40 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.
1import pygame
2import spring_mass as spm
3
4
5class ImageDrawer:
6 def __init__(self, image):
7 self.image = image
8
9 def __call__(self, screen, center, radius):
10 screen.blit(self.image, (center[0] - self.image.get_width() / 2,
11 center[1] - self.image.get_height() / 2))
12
48
49 image = pygame.image.load("../../assets/player/p1_walk03.png").convert_alpha()
50 image_drawer = ImageDrawer(image)
51 p5 = spm.PointMass((200, 300), (5, -5), world, radius=image.get_height() / 2,
52 mass=50, restitution=0.6, drawer=image_drawer)
53 actor_list.append(p5)
54
55 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.
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.
1import pygame
2import spring_mass as spm
3
4
5def create_point_mass(pos, world, fixed=False):
6 vel = (0, 0)
7 radius = 10
8 mass = 0.1
9 viscous = 0.01
10 restitution = 0.95
11 if fixed:
12 PointMassClass = spm.FixedPointMass
13 else:
14 PointMassClass = spm.PointMass
15 return PointMassClass(pos, vel, world, radius, mass, viscous, restitution,
16 spm.CircleDrawer("green", width=0))
17
18
19def create_spring(p1, p2, world):
20 spring_const = 0.02
21 natural_len = 5
22 return spm.Spring(p1, p2, world, spring_const, natural_len,
23 spm.LineDrawer("white", width=2))
24
25
26def main():
27 pygame.init()
28 width, height = 600, 400
29 screen = pygame.display.set_mode((width, height))
30 clock = pygame.time.Clock()
31
32 world = spm.World((width, height), dt=1.0, gravity_acc=(0, 0.5))
33 actor_list = []
34
35 p = []
36 sp = []
37 spacing = 30
38 p.append(create_point_mass((width/2, 0), world, fixed=True))
39 for k in range(1, 15):
40 p.append(create_point_mass((width/2, spacing * k), world))
41 sp.append(create_spring(p[k - 1], p[k], world))
42 actor_list += p + sp
43
44 while True:
45 frames_per_second = 60
46 clock.tick(frames_per_second)
47
48 should_quit = False
49 for event in pygame.event.get():
50 if event.type == pygame.QUIT:
51 should_quit = True
52 elif event.type == pygame.KEYDOWN:
53 if event.key == pygame.K_ESCAPE:
54 should_quit = True
55 elif event.key == pygame.K_SPACE:
56 sp[0].is_alive = False
57 if should_quit:
58 break
59
60 for a in actor_list:
61 a.update()
62 actor_list[:] = [a for a in actor_list if a.is_alive]
63
64 screen.fill(pygame.Color("black"))
65 for a in actor_list:
66 a.draw(screen)
67 pygame.display.update()
68
69 pygame.quit()
70
71
72if __name__ == "__main__":
73 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?
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.
1import pygame
2import spring_mass as spm
3
4
5class AppMain:
6 def __init__(self):
7 pygame.init()
8 width, height = 600, 400
9 self.screen = pygame.display.set_mode((width, height))
10
11 self.world = spm.World((width, height), dt=1.0, gravity_acc=(0, 0.1))
12 self.actor_list = []
13
14 def update(self):
15 for a in self.actor_list:
16 a.update()
17 self.actor_list[:] = [a for a in self.actor_list if a.is_alive]
18
19 def draw(self):
20 self.screen.fill(pygame.Color("black"))
21 for a in self.actor_list:
22 a.draw(self.screen)
23 pygame.display.update()
24
25 def run(self):
26 clock = pygame.time.Clock()
27
28 while True:
29 frames_per_second = 60
30 clock.tick(frames_per_second)
31
32 should_quit = False
33 for event in pygame.event.get():
34 if event.type == pygame.QUIT:
35 should_quit = True
36 elif event.type == pygame.KEYDOWN and event.key == pygame.K_ESCAPE:
37 should_quit = True
38 if should_quit:
39 break
40
41 self.update()
42 self.draw()
43
44 pygame.quit()
45
46
47if __name__ == "__main__":
48 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.
1import random
2import pygame
3import spring_mass as spm
4
5
6class ActorFactory:
7 def __init__(self, world, actor_list):
8 self.world = world
9 self.actor_list = actor_list
10
11 def create_point_mass(self, pos, fixed=False):
12 vel = (random.uniform(-10, 10), random.uniform(-10, 0))
13 mass = 10
14 radius = 10
15 viscous = 0.01
16 restitution = 0.95
17 if fixed:
18 PointMassClass = spm.FixedPointMass
19 color = "gray"
20 else:
21 PointMassClass = spm.PointMass
22 color = "green"
23 return PointMassClass(pos, vel, self.world, radius, mass, viscous,
24 restitution, spm.CircleDrawer(color, width=0))
25
26 def create_collision_resolver(self):
27 return spm.CollisionResolver(self.world, self.actor_list)
28
29 def create_boundary(self, name):
30 width, height = self.world.size
31 geometry = {"top": ((0, -1), (0, 0)),
32 "bottom": ((0, 1), (0, height)),
33 "left": ((-1, 0), (0, 0)),
34 "right": ((1, 0), (width, 0))}
35 normal, point_included = geometry[name]
36 return spm.Boundary(normal, point_included, self.world, self.actor_list)
37
38
39class AppMain:
40 def __init__(self):
41 pygame.init()
42 width, height = 600, 400
43 self.screen = pygame.display.set_mode((width, height))
44
45 self.world = spm.World((width, height), dt=1.0, gravity_acc=(0, 0.1))
46 self.actor_list = []
47 self.factory = ActorFactory(self.world, self.actor_list)
48
49 self.actor_list.append(self.factory.create_collision_resolver())
50 self.actor_list.append(self.factory.create_boundary("top"))
51 self.actor_list.append(self.factory.create_boundary("bottom"))
52 self.actor_list.append(self.factory.create_boundary("left"))
53 self.actor_list.append(self.factory.create_boundary("right"))
54
55 def update(self):
73 should_quit = False
74 for event in pygame.event.get():
75 if event.type == pygame.QUIT:
76 should_quit = True
77 elif event.type == pygame.KEYDOWN and event.key == pygame.K_ESCAPE:
78 should_quit = True
79 elif event.type == pygame.MOUSEBUTTONDOWN:
80 self.actor_list.append(self.factory.create_point_mass(event.pos))
81 if should_quit:
82 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.
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.
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.
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.
25
26 def create_spring(self, p1, p2):
27 spring_const = 0.01
28 natural_len = 20
29 return spm.Spring(p1, p2, self.world, spring_const, natural_len,
30 spm.LineDrawer("white", width=2))
31
32 def create_collision_resolver(self):
45class AppMain:
46 def __init__(self):
47 pygame.init()
48 width, height = 600, 400
49 self.screen = pygame.display.set_mode((width, height))
50
51 self.world = spm.World((width, height), dt=1.0, gravity_acc=(0, 0.1))
52 self.actor_list = []
53 self.factory = ActorFactory(self.world, self.actor_list)
54
55 self.actor_list.append(self.factory.create_collision_resolver())
56 self.actor_list.append(self.factory.create_boundary("top"))
57 self.actor_list.append(self.factory.create_boundary("bottom"))
58 self.actor_list.append(self.factory.create_boundary("left"))
59 self.actor_list.append(self.factory.create_boundary("right"))
60
61 self.point_mass_prev = None
62
63 def add_connected_point_mass(self, pos, button):
64 if button == 1:
65 fixed = False
66 elif button == 3:
67 fixed = True
68 else:
69 return
70 p = self.factory.create_point_mass(pos, fixed)
71 self.actor_list.append(p)
72
73 if self.point_mass_prev is not None:
74 sp = self.factory.create_spring(p, self.point_mass_prev)
75 self.actor_list.append(sp)
76 if pygame.key.get_pressed()[pygame.K_SPACE]:
77 self.point_mass_prev = p
78
79 def update(self):
97 should_quit = False
98 for event in pygame.event.get():
99 if event.type == pygame.QUIT:
100 should_quit = True
101 elif event.type == pygame.KEYDOWN and event.key == pygame.K_ESCAPE:
102 should_quit = True
103 elif event.type == pygame.KEYUP and event.key == pygame.K_SPACE:
104 self.point_mass_prev = None
105 elif event.type == pygame.MOUSEBUTTONDOWN:
106 self.add_connected_point_mass(event.pos, event.button)
107 if should_quit:
108 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.
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.
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.
1import pygame
2
3
4PgVector = pygame.math.Vector2
5
6
7class World:
8 def __init__(self, size, dt, gravity_acc):
9 self.size = size
10 self.dt = dt
11 self.gravity_acc = PgVector(gravity_acc)
12
13
14class CircleDrawer:
15 def __init__(self, color, width):
16 self.color = pygame.Color(color)
17 self.width = width
18
19 def __call__(self, screen, center, radius):
20 pygame.draw.circle(screen, self.color, center, radius, self.width)
21
22
23class LineDrawer:
24 def __init__(self, color, width):
25 self.color = pygame.Color(color)
26 self.width = width
27
28 def __call__(self, screen, pos1, pos2):
29 pygame.draw.line(screen, self.color, pos1, pos2, self.width)
30
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.
31
32def compute_gravity_force(mass, gravity_acc):
33 return mass * gravity_acc
34
35
36def compute_viscous_damping_force(viscous_damping, vel):
37 return -viscous_damping * vel
38
39
40def integrate_symplectic(pos, vel, force, mass, dt):
41 vel_new = vel + force / mass * dt
42 pos_new = pos + vel_new * dt
43 return pos_new, vel_new
44
45
46class PointMass:
47 def __init__(self, pos, vel, world, radius=10.0, mass=10.0,
48 viscous_damping=0.01, restitution=0.95, drawer=None):
49 self.is_alive = True
50 self.world = world
51 self.drawer = drawer or CircleDrawer("blue", 1)
52
53 self.pos = PgVector(pos)
54 self.vel = PgVector(vel)
55 self.radius = radius
56 self.mass = mass
57 self.viscous_damping = viscous_damping
58 self.restitution = restitution
59
60 self.total_force = PgVector((0, 0))
61
62 def update(self):
63 self.generate_force()
64 self.move()
65 self.total_force = PgVector((0, 0))
66
67 def draw(self, screen):
68 self.drawer(screen, self.pos, self.radius)
69
70 def receive_force(self, force):
71 self.total_force += PgVector(force)
72
73 def generate_force(self):
74 force_g = compute_gravity_force(self.mass, self.world.gravity_acc)
75 force_v = compute_viscous_damping_force(self.viscous_damping, self.vel)
76 self.receive_force(force_g + force_v)
77
78 def move(self):
79 self.pos, self.vel = \
80 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).
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.
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.
45class AppMain:
46 def __init__(self):
47 pygame.init()
48 width, height = 600, 400
49 self.screen = pygame.display.set_mode((width, height))
50
51 self.world = spm.World((width, height), dt=1.0, gravity_acc=(0, 0.1))
52 self.actor_list = []
53 self.factory = ActorFactory(self.world, self.actor_list)
54
55 # self.actor_list.append(self.factory.create_collision_resolver())
56 # self.actor_list.append(self.factory.create_boundary("top"))
57 # self.actor_list.append(self.factory.create_boundary("bottom"))
58 # self.actor_list.append(self.factory.create_boundary("left"))
59 # self.actor_list.append(self.factory.create_boundary("right"))
60
61 self.point_mass_prev = None
62
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.
82
83class FixedPointMass(PointMass):
84 def __init__(self, pos, vel, world, radius=10, mass=10,
85 viscous_damping=0.01, restitution=0.95, drawer=None):
86 super().__init__(pos, vel, world, radius, mass,
87 viscous_damping, restitution, drawer)
88 self.vel, self.mass = PgVector((0, 0)), 1e9
89
90 def move(self):
91 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.
93
94def compute_restoring_force(pos1, pos2, spring_const, natural_len):
95 if pos1 == pos2:
96 return None
97 vector12 = pos2 - pos1
98 distance = vector12.magnitude()
99 unit_vector12 = vector12 / distance
100 f1 = unit_vector12 * spring_const * (distance - natural_len)
101 return f1
102
103
104class Spring:
105 def __init__(self, point_mass1, point_mass2, world,
106 spring_const=0.01, natural_len=0.0, drawer=None):
107 self.is_alive = True
108 self.world = world
109 self.drawer = drawer or LineDrawer("blue", 1)
110
111 self.p1 = point_mass1
112 self.p2 = point_mass2
113 self.spring_const = spring_const
114 self.natural_len = natural_len
115
116 def update(self):
117 if not (self.p1.is_alive and self.p2.is_alive):
118 self.is_alive = False
119 return
120 self.generate_force()
121
122 def draw(self, screen):
123 self.drawer(screen, self.p1.pos, self.p2.pos)
124
125 def generate_force(self):
126 f1 = compute_restoring_force(self.p1.pos, self.p2.pos, self.spring_const, self.natural_len)
127 if f1 is None:
128 return
129 self.p1.receive_force(f1)
130 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:
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.
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.
1import pytest
2import spring_mass as spm
3from spring_mass import PgVector
4
5
6@pytest.mark.parametrize("pos1v, pos2v, k, nl, f1_expected, expects_none", [
7 ((0, 0), (10, 0), 2, 0, (2 * 10, 0), False),
8 ((10, 0), (0, 0), 2, 0, (-2 * 10, 0), False),
9 ((0, 0), (10, 0), 2, 7, (2 * (10 - 7), 0), False),
10 ((10, 0), (0, 0), 2, 7, (-2 * (10 - 7), 0), False),
11 ((0, 0), (10, 0), 2, 13, (-2 * (13 - 10), 0), False),
12 ((10, 0), (0, 0), 2, 13, (2 * (13 - 10), 0), False),
13
14 ((0, 0), (0, 10), 3, 0, (0, 3 * 10), False),
15 ((2, 10), (2, 5), 4, 0, (0, -5 * 4), False),
16 ((0, 0), (0, 10), 3, 5, (0, 3 * 5), False),
17 ((0, 0), (0, 10), 3, 12, (0, -3 * 2), False),
18
19 ((0, 0), (0, 0), 0, 0, (0, 0), True), # equal positions
20 ((12, 5), (12, 5), 0, 0, (0, 0), True), # equal positions
21])
22def test_restoring_force(pos1v, pos2v, k, nl, f1_expected, expects_none):
23 pos1 = PgVector(pos1v)
24 pos2 = PgVector(pos2v)
25 f1 = spm.compute_restoring_force(pos1, pos2, k, nl)
26 if expects_none:
27 assert f1 is None
28 else:
29 assert f1 == PgVector(f1_expected)
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.
132
133class FragileSpring(Spring):
134 def __init__(self, point_mass1, point_mass2, world,
135 spring_const=0.01, natural_len=0.0, drawer=None,
136 break_threshold=1e9):
137 super().__init__(point_mass1, point_mass2, world, spring_const,
138 natural_len, drawer)
139 self.break_threshold = break_threshold
140
141 def generate_force(self):
142 f1 = compute_restoring_force(self.p1.pos, self.p2.pos, self.spring_const, self.natural_len)
143 if f1 is None:
144 return
145 self.p1.receive_force(f1)
146 self.p2.receive_force(-f1)
147 if f1.magnitude() > self.break_threshold:
148 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.
26 def create_spring(self, p1, p2):
27 spring_const = 0.01
28 natural_len = 20
29 break_threshold = 5.0
30 return spm.FragileSpring(p1, p2, self.world, spring_const, natural_len,
31 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.
150
151def is_point_mass(actor):
152 return isinstance(actor, PointMass)
153
154
155def compute_impact_force_between_points(p1, p2, dt):
156 if (p1.pos - p2.pos).magnitude() > p1.radius + p2.radius:
157 return None
158 if p1.pos == p2.pos:
159 return None
160 normal = (p2.pos - p1.pos).normalize()
161 v1 = p1.vel.dot(normal)
162 v2 = p2.vel.dot(normal)
163 if v1 < v2:
164 return None
165 e = p1.restitution * p2.restitution
166 m1, m2 = p1.mass, p2.mass
167 f1 = normal * (-(e + 1) * v1 + (e + 1) * v2) / (1/m1 + 1/m2) / dt
168 return f1
169
170
171class CollisionResolver:
172 def __init__(self, world, actor_list, target_condition=None, drawer=None):
173 self.is_alive = True
174 self.world = world
175 self.drawer = drawer or (lambda surface: None)
176
177 self.actor_list = actor_list
178 if target_condition is None:
179 self.target_condition = is_point_mass
180 else:
181 self.target_condition = target_condition
182
183 def update(self):
184 self.generate_force()
185
186 def draw(self, surface):
187 self.drawer(surface)
188
189
190 def generate_force(self):
191 plist = [a for a in self.actor_list if self.target_condition(a)]
192 n = len(plist)
193 for i in range(n):
194 for j in range(i + 1, n):
195 p1, p2 = plist[i], plist[j]
196 f1 = compute_impact_force_between_points(p1, p2, self.world.dt)
197 if f1 is None:
198 continue
199 p1.receive_force(f1)
200 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
and with the restitution coefficient e,
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
By solving these equations, we get the following force:
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.
54 self.factory = ActorFactory(self.world, self.actor_list)
55
56 self.actor_list.append(self.factory.create_collision_resolver())
57 # self.actor_list.append(self.factory.create_boundary("top"))
58 # self.actor_list.append(self.factory.create_boundary("bottom"))
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.
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.
31
32@pytest.mark.parametrize("p1v, p2v, expects_none", [
33 (((0, 0), (5, 0), 20, 10, 1.0), # (pos, vel, r, m, e)
34 ((20, 0), (-5, 0), 20, 10, 1.0), False),
35 (((0, 0), (5, 0), 20, 15, 0.2),
36 ((20, 0), (-5, 0), 20, 30, 0.8), False),
37
38 (((0, 0), (5, 0), 20, 10, 1.0),
39 ((120, 0), (-5, 0), 20, 10, 1.0), True), # not in contact
40 (((0, 0), (5, 0), 20, 10, 1.0),
41 ((0, 0), (-5, 0), 20, 10, 1.0), True), # equal positions
42 (((0, 0), (-5, 0), 20, 10, 1.0),
43 ((20, 0), (5, 0), 20, 10, 1.0), True), # moving apart
44])
45def test_impact_force_between_points(p1v, p2v, expects_none):
46 dt, g = 1.0, 0
47 w = spm.World((600, 400), dt, g)
48 pos1, vel1, r1, m1, e1 = p1v
49 pos2, vel2, r2, m2, e2 = p2v
50
51 p1 = spm.PointMass(pos1, vel1, w,
52 radius=r1, mass=m1, restitution=e1)
53 p2 = spm.PointMass(pos2, vel2, w,
54 radius=r2, mass=m2, restitution=e2)
55 f1 = spm.compute_impact_force_between_points(p1, p2, dt)
56 if expects_none:
57 assert f1 is None
58 return
59
60 vel1_old = PgVector(vel1)
61 vel2_old = PgVector(vel2)
62 vel1_new = vel1_old + f1 / m1 * dt
63 vel2_new = vel2_old - f1 / m2 * dt
64 e = e1 * e2
65 assert vel1_new - vel2_new == -e * (vel1_old - vel2_old)
66 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.
202
203def compute_impact_force_by_fixture(p, normal, point_included, dt):
204 invasion = normal.dot(p.pos - point_included)
205 if invasion + p.radius > 0 and normal.dot(p.vel) > 0:
206 e = p.restitution
207 v = normal.dot(p.vel)
208 m = p.mass
209 f = normal * (-(e + 1) * v) * m / dt
210 else:
211 f = None
212 return f
213
214
215class Boundary:
216 def __init__(self, normal, point_included, world, actor_list,
217 target_condition=None, drawer=None):
218 self.is_alive = True
219 self.world = world
220 self.drawer = drawer or (lambda surface: None)
221
222 self.normal = PgVector(normal).normalize()
223 self.point_included = PgVector(point_included)
224 self.actor_list = actor_list
225 if target_condition is None:
226 self.target_condition = is_point_mass
227 else:
228 self.target_condition = target_condition
229
230 def update(self):
231 self.generate_force()
232
233 def draw(self, surface):
234 self.drawer(surface)
235
236
237 def generate_force(self):
238 plist = [a for a in self.actor_list if self.target_condition(a)]
239 for p in plist:
240 f = compute_impact_force_by_fixture(p, self.normal, self.point_included, self.world.dt)
241 if f is None:
242 continue
243 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.
56 self.actor_list.append(self.factory.create_collision_resolver())
57 self.actor_list.append(self.factory.create_boundary("top"))
58 self.actor_list.append(self.factory.create_boundary("bottom"))
59 self.actor_list.append(self.factory.create_boundary("left"))
60 self.actor_list.append(self.factory.create_boundary("right"))
61
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.
68
69@pytest.mark.parametrize("pv, n, pinc, vel_expected, expects_none", [
70 (((600, 200), (5, 3), 20, 10, 1.0), # (pos, vel, r, m, e)
71 (1, 0), (600, 0), (-5, 3), False),
72 (((15, 200), (-5, 3), 20, 10, 0.8),
73 (-1, 0), (0, 0), (0.8 * 5, 3), False),
74 (((300, 400), (-5, 3), 20, 10, 0.6),
75 (0, 1), (0, 400), (-5, - 0.6 * 3), False),
76
77 (((600, 200), (-5, 3), 20, 10, 1.0),
78 (1, 0), (600, 0), (0, 0), True), # not moving toward warll
79 (((580, 200), (5, 3), 20, 10, 1.0),
80 (1, 0), (600, 0), (0, 0), True), # not in contact
81])
82def test_impact_force_by_fixture(pv, n, pinc, vel_expected, expects_none):
83 dt, g = 1.0, 0
84 w = spm.World((600, 400), dt, g)
85 pos, vel, r, m, e = pv
86 p = spm.PointMass(pos, vel, w, radius=r, mass=m, restitution=e)
87 normal = PgVector(n)
88 point_included = PgVector(pinc)
89 f = spm.compute_impact_force_by_fixture(p, normal, point_included, dt)
90 if expects_none:
91 assert f is None
92 return
93
94 vel_old = PgVector(vel)
95 vel_new = vel_old + f / m * dt
96 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.
1import pygame
2
3
4PgVector = pygame.math.Vector2
5
6
7class World:
8 def __init__(self, size, dt, gravity_acc):
9 self.size = size
10 self.dt = dt
11 self.gravity_acc = PgVector(gravity_acc)
12
13
14class CircleDrawer:
15 def __init__(self, color, width):
16 self.color = pygame.Color(color)
17 self.width = width
18
19 def __call__(self, screen, center, radius):
20 pygame.draw.circle(screen, self.color, center, radius, self.width)
21
22
23class LineDrawer:
24 def __init__(self, color, width):
25 self.color = pygame.Color(color)
26 self.width = width
27
28 def __call__(self, screen, pos1, pos2):
29 pygame.draw.line(screen, self.color, pos1, pos2, self.width)
30
31
32def compute_gravity_force(mass, gravity_acc):
33 return mass * gravity_acc
34
35
36def compute_viscous_damping_force(viscous_damping, vel):
37 return -viscous_damping * vel
38
39
40def integrate_symplectic(pos, vel, force, mass, dt):
41 vel_new = vel + force / mass * dt
42 pos_new = pos + vel_new * dt
43 return pos_new, vel_new
44
45
46class PointMass:
47 def __init__(self, pos, vel, world, radius=10.0, mass=10.0,
48 viscous_damping=0.01, restitution=0.95, drawer=None):
49 self.is_alive = True
50 self.world = world
51 self.drawer = drawer or CircleDrawer("blue", 1)
52
53 self.pos = PgVector(pos)
54 self.vel = PgVector(vel)
55 self.radius = radius
56 self.mass = mass
57 self.viscous_damping = viscous_damping
58 self.restitution = restitution
59
60 self.total_force = PgVector((0, 0))
61 self.message_list = []
62
63 def update(self):
64 self.generate_force()
65 self.move()
66 self.total_force = PgVector((0, 0))
67 self.message_list.clear()
68
69 def draw(self, screen):
70 self.drawer(screen, self.pos, self.radius)
71
72 def receive_force(self, force):
73 self.total_force += PgVector(force)
74
75 def receive_message(self, msg):
76 self.message_list.append(msg)
77
78 def generate_force(self):
79 force_g = compute_gravity_force(self.mass, self.world.gravity_acc)
80 force_v = compute_viscous_damping_force(self.viscous_damping, self.vel)
81 self.receive_force(force_g + force_v)
82
83 def move(self):
84 self.pos, self.vel = \
85 integrate_symplectic(self.pos, self.vel, self.total_force, self.mass, self.world.dt)
86
87 for msg in self.message_list:
88 if msg["type"] == "floor_hit" and self.vel.y > 0:
89 # constrain y on or above floor
90 self.pos.y = msg["y"] - self.radius
91
92
93class FixedPointMass(PointMass):
94 def __init__(self, pos, vel, world, radius=10.0, mass=10.0,
95 viscous_damping=0.01, restitution=0.95, drawer=None):
96 super().__init__(pos, vel, world, radius, mass,
97 viscous_damping, restitution, drawer)
98 self.vel, self.mass = PgVector((0, 0)), 1e9
99
100 def move(self):
101 pass
102
103
104def compute_restoring_force(pos1, pos2, spring_const, natural_len):
105 if pos1 == pos2:
106 return None
107 vector12 = pos2 - pos1
108 distance = vector12.magnitude()
109 unit_vector12 = vector12 / distance
110 f1 = unit_vector12 * spring_const * (distance - natural_len)
111 return f1
112
113
114class Spring:
115 def __init__(self, point_mass1, point_mass2, world,
116 spring_const=0.01, natural_len=0.0, drawer=None):
117 self.is_alive = True
118 self.world = world
119 self.drawer = drawer or LineDrawer("blue", 1)
120
121 self.p1 = point_mass1
122 self.p2 = point_mass2
123 self.spring_const = spring_const
124 self.natural_len = natural_len
125
126 def update(self):
127 if not (self.p1.is_alive and self.p2.is_alive):
128 self.is_alive = False
129 return
130 self.generate_force()
131
132 def draw(self, screen):
133 self.drawer(screen, self.p1.pos, self.p2.pos)
134
135 def generate_force(self):
136 f1 = compute_restoring_force(self.p1.pos, self.p2.pos, self.spring_const, self.natural_len)
137 if f1 is None:
138 return
139 self.p1.receive_force(f1)
140 self.p2.receive_force(-f1)
141
142
143class FragileSpring(Spring):
144 def __init__(self, point_mass1, point_mass2, world,
145 spring_const=0.01, natural_len=0.0, drawer=None,
146 break_threshold=1e9):
147 super().__init__(point_mass1, point_mass2, world, spring_const,
148 natural_len, drawer)
149 self.break_threshold = break_threshold
150
151 def generate_force(self):
152 f1 = compute_restoring_force(self.p1.pos, self.p2.pos, self.spring_const, self.natural_len)
153 if f1 is None:
154 return
155 self.p1.receive_force(f1)
156 self.p2.receive_force(-f1)
157 if f1.magnitude() > self.break_threshold:
158 self.is_alive = False
159
160
161def is_point_mass(actor):
162 return isinstance(actor, PointMass)
163
164
165def compute_impact_force_between_points(p1, p2, dt):
166 if (p1.pos - p2.pos).magnitude() > p1.radius + p2.radius:
167 return None
168 if p1.pos == p2.pos:
169 return None
170 normal = (p2.pos - p1.pos).normalize()
171 v1 = p1.vel.dot(normal)
172 v2 = p2.vel.dot(normal)
173 if v1 < v2:
174 return None
175 e = p1.restitution * p2.restitution
176 m1, m2 = p1.mass, p2.mass
177 f1 = normal * (-(e + 1) * v1 + (e + 1) * v2) / (1/m1 + 1/m2) / dt
178 return f1
179
180
181class CollisionResolver:
182 def __init__(self, world, actor_list, target_condition=None, drawer=None):
183 self.is_alive = True
184 self.world = world
185 self.drawer = drawer or (lambda surface: None)
186
187 self.actor_list = actor_list
188 if target_condition is None:
189 self.target_condition = is_point_mass
190 else:
191 self.target_condition = target_condition
192
193 def update(self):
194 self.generate_force()
195
196 def draw(self, surface):
197 self.drawer(surface)
198
199
200 def generate_force(self):
201 plist = [a for a in self.actor_list if self.target_condition(a)]
202 n = len(plist)
203 for i in range(n):
204 for j in range(i + 1, n):
205 p1, p2 = plist[i], plist[j]
206 f1 = compute_impact_force_between_points(p1, p2, self.world.dt)
207 if f1 is None:
208 continue
209 p1.receive_force(f1)
210 p2.receive_force(-f1)
211
212
213def compute_impact_force_by_fixture(p, normal, point_included, dt):
214 invasion = normal.dot(p.pos - point_included)
215 if invasion + p.radius > 0 and normal.dot(p.vel) > 0:
216 e = p.restitution
217 v = normal.dot(p.vel)
218 m = p.mass
219 f = normal * (-(e + 1) * v) * m / dt
220 else:
221 f = None
222 return f
223
224
225class Boundary:
226 def __init__(self, normal, point_included, world, actor_list,
227 target_condition=None, drawer=None):
228 self.is_alive = True
229 self.world = world
230 self.drawer = drawer or (lambda surface: None)
231
232 self.normal = PgVector(normal).normalize()
233 self.point_included = PgVector(point_included)
234 self.actor_list = actor_list
235 if target_condition is None:
236 self.target_condition = is_point_mass
237 else:
238 self.target_condition = target_condition
239
240 def update(self):
241 self.generate_force()
242
243 def draw(self, surface):
244 self.drawer(surface)
245
246
247 def is_floor(self):
248 return self.normal == PgVector((0, 1))
249
250 def generate_force(self):
251 plist = [a for a in self.actor_list if self.target_condition(a)]
252 for p in plist:
253 f = compute_impact_force_by_fixture(p, self.normal, self.point_included, self.world.dt)
254 if f is None:
255 continue
256 p.receive_force(f)
257 if self.is_floor():
258 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.
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.
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.