7. バネ・質点系の力学シミュレーション¶
前章までの運動学シミュレーションでは,質点どうしの間には相互作用がありませんでした.この章では,力学シミュレーションを題材として,オブジェクト間にもう少し複雑な相互作用があるプログラムの実際を学びます.
難しそうに聞こえるかもしれませんが,これまで学んできたクラスの継承やオブジェクトの合成などの考え方を活用することで,単純な要素の組み合わせで複雑なシステムを構築できることを見ていきます.
- 左クリック: 質点を生成
- 右クリック: 固定質点を生成
- スペースキーを押しながらクリック: バネでつながれた質点を生成
- ESCキー: 終了
7.1. まず試してみる¶
この章で扱う力学シミュレーションは,質点,バネ,質点間の衝突,質点と壁の衝突を扱います.これらは spring_mass モジュールとして利用できるようにします.
spring_mass モジュールを作る話は後回しにして,まずはユーザとして使ってみるところから始めます.
これまで同様 projects
フォルダの中に spring_mass
というフォルダを作成し,VSCode でそのフォルダを開いた状態にしてください.そこに spring_mass.py
というファイルを作成して,この章の最後の spring_mass.py
(ver 18.0) のコードをコピーペーストして保存しておいてください.この時点で VSCode のタイトルバーには
spring_mass.py - spring_mass - VSCodium (Computer Seminar I)
と表示されているはずです.
次に,spring_mass モジュールを使う側のプログラムとして
simple_main.py
ファイルを作成し,以下の内容を入力してください.
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()
基本構造は particle.py
の ver 1.0 と一緒です.
2行目で spring_mass モジュールを import
しています.名前が長いので
spm
という名前で使えるようにしています.
以降は main
関数だけのプログラムで,いつも通りの空っぽの雛形に,
World
型オブジェクトと actor_list
というリストを用意しています
(11-12行目).
World
の役割はこれまでと同じですが,初期化の際の引数の仕様が少し変わっています.第1引数には画面の幅と高さのタプルを渡します.第2引数に単位時間 dt
,第3引数は重力加速度ベクトルをタプルで渡します (これまでは y
成分だけのスカラ値でした).
この章では,質点やバネなどの構成要素を総称して Actor と呼ぶことにします.これまでの particle_list に相当するものを一般化して actor_list
という名前にしています.
28-30行目で actor_list
内の Actor が update
され,is_alive
でないものは削除されます.33-34行目で各 Actor が draw
されます.
真っ黒な画面が出て,Esc キーで終了できれば OK です.
Source Control サイドバーで git リポジトリを初期化し,simple_main.py
をコミットしておいてください.
spring_mass.py
の方も一緒にコミットしてもよいですし,しなくてもよいです.
7.1.1. 質点を配置¶
質点を1個生成します.前章までの Particle
は属性として位置,速度,半径
(と色) を持っていましたが,この章では力学的な相互作用を実現するため,質量,粘性抵抗係数,反発係数を持たせます.これらは初期化時に引数として渡すことで設定します.力学特性を持つという点で前章までとは大きく変わるので,クラス名も新たに PointMass
としました.
World
型オブジェクトへの参照を持つのもこれまでと同じです.一つ大きく変わるのは,描画処理を Strategy パターンで入れ替え可能になっている点です.spring_mass モジュールに CircleDrawer
という円を描画するオブジェクトのクラスが定義されています.色と線幅 (0 なら塗りつぶし)
を指定して生成します.
描画処理を入れ替え可能とすることで,単に円を描くだけでなく,もっと凝ったグラフィクスを描画したり,画像を貼り付けたりすることができるようになります.具体例は後で見ることにします.
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:
18行目では緑色,塗りつぶしの circle_drawer
を生成し,19行目の
PointMass
型オブジェクトの生成時に渡しています.生成された PointMass
型オブジェクトは actor_list
に append
されています.
なお,radius
以降の引数にはデフォルト値が設定されています.具体的な値については spring_mass.py
内の定義を見てください.
7.1.2. 壁・床を配置¶
particle.py
では壁で跳ね返るのは Particle
自身の責務 (正確には
Particle
が属性として保持しているストラテジオブジェクトの責務) でしたが,この章では,質点に力を及ぼす壁 (Boundary
型のオブジェクト) を生成し,質点は力を受けて結果的に跳ね返ることにします.
Boundary
型オブジェクトを生成する際は,第1引数に法線ベクトル,第2引数に壁面上に含まれる任意の1点の座標,第3引数に World
型オブジェクト,第4
引数に衝突対象になる質点を含むリストを渡します (ここに質点以外が含まれる場合は無視されます).法線ベクトルは,壁の中に進む向きとします.
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:
22行目で作られるのは (width, 0) を通り右向き法線を持つ境界なので,画面右端の壁になります.23行目で床を作ります.左側の壁や天井なども同様に作れます.
これらを actor_list
に入れておけば,対応する境界での跳ね返りが発生します.actor_list
から外せば跳ね返らず外に飛び去っていきます.いろいろ試してみてください.
これを解消するには,「壁の中」に一定距離以上入り込んだときは力を発生しないなどの対処が考えられます.ただしうまく調整しないと質点の速度によっては「すり抜け」が発生するかも知れません.
また,ここで作った壁には「表」と「裏」があることに注意してください.表からも裏からも衝突できる壁が欲しい場合は,表向きの壁と裏向きの壁 (それぞれ法線方向が互いと逆) の 2 枚を同じところに置くなどの対処が必要です.
7.1.3. バネを配置¶
バネはクラス Spring
で表します.初期化時に,両端の質点 2 つと
World
型オブジェクトのほか,バネ定数,自然長,描画ストラテジオブジェクトを渡します.
直線を描画するオブジェクトのクラスとして LineDrawer
が用意されています.
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:
ここでは適当な 2 つの質点を生成し,バネでつないでいます.重力があるとバネ自体の動作がわかりにくいかも知れません.11行目の重力加速度を適宜
(0, 0)
にするなどしてみてください.
質点間の衝突はまだ発生しない点に注意してください.
PointMass
を継承した FixedPointMass
クラスも用意されています.固定端につながれたバネ質量は以下のように配置できます.
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. 質点間の衝突解決¶
質点間の衝突を発生させるには,CollisionResolver
型のオブジェクトを
actor_list
に入れます.
初期化の際は,World
型オブジェクトと,衝突判定の対象を含むリストを渡します.ここに含まれる質点の相互間で衝突が発生します.
37
38 actor_list.append(spm.CollisionResolver(world, actor_list))
39
40 while True:
7.1.5. 描画処理のカスタマイズ¶
描画処理クラスを自作する例を示します.まずは spring_mass モジュールが用意している CircleDrawer
と LineDrawer
の定義を見るのが早いです:
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)
つまり,PointMass
生成時に渡す描画オブジェクトは,引数 screen
, center
, radius
を受け取る呼び出し可能オブジェクトであり,
Spring
生成時に渡す方は screen
と両端の位置 pos1
, pos2
を受け取る呼び出し可能オブジェクトです.
これに倣って CircleDrawer
の代わりに適当な画像を表示する
ImageDrawer
を作ってみましょう.
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:
ImageDrawer
の定義は 5-11行目の通りです.初期化時に画像
(Surface
型オブジェクト) を受け取り,__call__
される度に blit
で貼り付けます.
49-50行目では,hello_pygame.py
で使ったプレイヤ画像を読み込み,ImageDrawer
型オブジェクトを生成しています.後は CircleDrawer
と同じ使い方です.
なお,画像の形がどんなものでも,衝突計算は円形として行われます.
このテキストでは,円どうしの衝突と,円と直線の衝突のみを扱っています.これらの判定は簡単だからです.もっと複雑な形状の衝突を考えようとすると,一般論はとても難しく,それだけで本が一冊書けてしまうくらいです.
円と直線以外で比較的扱いやすいものとしては,縦横が座標軸方向と揃っている軸平行長方形 (axis-aligned rectangle, aa-rectangle) や,長方形の両端に半円をくっつけたカプセル (capsule) などがあります.ゲームなどのように厳密性が求められない場合は,複雑な形状もこれらで近似してしまうことも多いです.興味のある人は調べてみてください.
なお,pygame には軸平行長方形の衝突判定を行う sprite というモジュールが用意されています.
7.2. Slinky Drop シミュレーション¶
少し応用っぽいものを作ってみます.Slinky Drop という理科実験の話を聞いたことはあるでしょうか.以下の動画で紹介されています.
この動画は実験結果が判明する直前で終わります.結果が気になる人は,続きの動画を探すか,またはせっかくなので自分でシミュレーションしてみてはどうでしょう.
simple_main.py
の最初のバージョン (ver 1.0) をコピーして,
slinky_drop_main.py
を作り,バネと質量を追加してスリンキーを模擬します.
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()
- 5-16行目,19-23行目:
- 質点の生成,バネの生成を行う関数を用意しておきます.メインプログラム側では,これらを呼び出すことで細かい引数の設定などに煩わされることなくオブジェクトを生成できます.このようにオブジェクト生成を行う関数をファクトリ関数と呼ぶことがあります.
- 35-42行目
- 画面中央の一番上に固定質点を置き,そこから15個の質点をバネで直列につなぎます.
- 52-56行目
- Space キーを押すと,最上部のバネの
is_alive
属性がFalse
になり,手を離したことを模擬します.
さて,予想通りの動きだったでしょうか?
sp[0].is_alive = False
としているところは「属性はクラスの外から書き換えるのは避ける」というこのテキストの方針に違反してませんか?真面目に方針を守るなら,Spring
クラスに自身の is_alive
属性を False
にするメソッド (例えば名前は disappear
とか)
を設けてそれを呼ぶべきです.あるいは Actor はすべてこのメソッドを持つべきかも知れません.
「is_alive
を False
にするだけの処理をわざわざメソッドにする必要があるのか?」と考える人も多いと思います.実際,Python プログラマの場合はこのくらいは気にせずにクラス外から直接書き換える人も多く,その方が Python らしいという人もいます.
このテキストで「属性を外から書き換えない」方針をできるだけ守ってきたのは,オブジェクトに責務を任せる (オブジェクトの状態を変更するのはオブジェクト自身にやらせる) という考え方に慣れてもらうためです.その考え方が十分に理解できるようになったなら,実際のコードの書き方のルールは自分で(あるいはチームで) 決めればよいと思います.
7.3. インタラクティブ放物シミュレーション¶
もう少し複雑な例として,この章の冒頭にアニメーションを示した,マウス・キー操作でインタラクティブに質量やバネを配置できるものを作ります.
まずは空っぽのフレームワークです.ちょっと複雑になるのでメインプログラム側もクラスとして作ります.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. 質点と壁と衝突処理¶
まず,質点と壁を生成し,衝突処理を発動させる部分だけ作ります.これらはオブジェクトを配置するだけなので,比較的簡単です.
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
- 6-36行目
利便性のため,ファクトリ関数群をメソッドとして持つクラス
ActorFactory
を作っておきます.まずこのオブジェクトにworld
やactor_list
を保持させておき,Actor を生成するときに内部で渡してやるようにすれば,メインプログラム側でいちいちこれらを引数として渡す必要がなくなります.29-36行目では
Boundary
のファクトリを定義していますが,上下左右の境界を法線ベクトルと点の座標で指定するのは煩雑なので,名前で指定できるようにしています."top"
,"bottom"
,"left"
,"right"
という文字列をキーとして,法線ベクトルと点の座標のタプルを返す辞書geometry
を作り,必要な値を入れておいてオブジェクト生成時に使います.- 47-53行目
- ファクトリオブジェクト
factory
を生成し,これを使って衝突処理オブジェクトと境界オブジェクトを生成し,actor_list
に加えます.まあまあすっきり書けるようになったと思います. - 79-80行目
- マウスボタンが押されたら質点を生成し,
actor_list
に加えます.この時点では,どのボタンかに関わらず非固定の質点が生成されます.
PointMassClass
という変数に「クラス名」を代入していますが,こんなことしていいんですか?過去の Q&A で何度か言及しましたが,Python ではクラスもオブジェクトの一種として扱えて,変数に割り当てたり,関数に引数として渡したりできます.23行目のように変数を介してインスタンスを生成することもできます.
ここでは PointMass
と FixedPointMass
を生成するときの引数は全く同じでよいので,if
と else
のそれぞれで生成するより,このように書く方が簡潔です.
AppMain
のメソッドにしてしまえば,そもそも world
や actor_list
を渡す必要なくなって楽なんじゃないですか?AppMain
クラスの肥大化を防ぐために,外部に出せるものはできるだけ外部に出すことにしています.何度も話しているように,プログラムの大規模化に立ち向かうときの鍵は,同時に考慮しなくてはならないものを減らすことです.クラスはできるだけ小さく保つ方が,同時に考慮すべきものが少なくなります.
ActorFactory
が持つメソッドは,AppMain
の属性を変更しないように作られています.こういうものはできるだけ外に追い出す方がよいです.
一方,AppMain
の属性を更新するような処理は迂闊にクラス外に追い出せません (このテキストでは,属性はできるだけクラス外からは変更しない方針を採ってきたのを思い出してください).
ActorFactory
が AppMain
の actor_list
属性への参照を保持しているからといって,「そうだ,create_point_mass
メソッドの中で actor_list
に append
してしまえば一石二鳥だ!」などとは考えない方がよいです.
よく知られている経験則として「オブジェクトを作るクラスとそのオブジェクトを使うクラスは分けておく方がよい」という設計指針が挙げられます.ファクトリという「オブジェクトを生成するだけ」の関数やクラスを作るのは,そういう指針の一環であるといえます.
少しだけメリットが現れているのは,world
や actor_list
のように共通で必要となるようなものを毎回渡さずに済む点です.
もっと有難みが出るのは,生成したいオブジェクトの組み合わせをごっそり入れ替えられるようにしたい場合です.次の節で
ActorFactory
にバネを生成するメソッドを追加します.バネや質点のパラメータ群はこの章では 1 種類に固定していますが,目的に応じて異なるバネマスセット (例えば重くて硬いセットとか,柔らかく自然長の長いセットとか) を用意したい場合もあり得ます.パラメータを外部ファイルから読み込みたい場合などもあるでしょう.
そのような場合,それぞれのバネマスセットを ActorFactory
の派生クラスとして定義しておけば,47行目で作る factory
を変えるだけで簡単に切り替えられるようになります.
7.3.2. バネを追加¶
あとはバネを作る方法を加えれば完成です.Space キーを押している間は,生成した質点がバネで順につながれるようにします.
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
- 26-30行目
- ファクトリに
Spring
を作るメソッドを加えておきます. - 61行目
- 質点を順にバネでつなぐためには,一つ前に作った質点を覚えておく必要があります.ちょっと安直ですが,
AppMain
の属性point_mass_prev
として覚えておくことにします.初期値はNone
(覚えてない) です. - 63-77行目
マウスボタンが押されたときの処理は,ちょっと複雑なのでメソッドに分けておくことにします.
まず64-71行目で,ボタン番号に応じて非固定または固定の質点を生成し,
actor_list
に入れておきます.次の73-75行目で,一つ前の質点を覚えているなら,さっき作った質点とバネでつなぎます.
最後に,Space キーが押されているなら,さっき作った質点を
point_mass_prev
属性として記憶しておきます.- 103-104行目
- 忘れてはいけないのは,Space キーから手を離したときに
point_mass_prev
をNone
にリセットしておくことです.これをしないと何が起きるか想像できるでしょうか.この2行をコメントアウトして答え合わせしてみてください. - 106行目
- マウスボタンが押されたときは
add_connected_point_mass
に処理をさせます.
run
と
add_connected_point_mass
に分散しているのがわかりにくいです.この点に限らず,このテキスト全体を通して,キー・マウス入力処理は抽象化が足りていません.あまり凝った構造にしても学びにくいと考えたためです.
本来は,キー・マウス入力に関する処理は AppMain
から独立させた上で,必要に応じてそちらから AppMain
のメソッドを呼び出させたり,AppMain
側からの状態問い合わせを受け付けたりするように設計する方がよいでしょう.
if
- elif
- elif
... を羅列するのもあまり読みやすいコードとはいえません.イベント処理を行うクラスをイベントの種類ごとに用意して,event.type
をキーとする辞書に入れておくとか,「Chain of Responsibility」と呼ばれるデザインパターンを利用するとか,改善方法はいろいろと考えられます.
7.4. spring_mass モジュールを作る¶
アプリケーション側の事例はここまでとして,これまで天下りで与えていた spring_mass モジュールを見ていきます.
7.4.1. 基本設計¶
まず基本的な構造として,Actor はメソッド update
と draw
を持ち,
is_alive
属性を持つものとします.
update
が呼ばれると,Actor は自身が力を及ぼすべき相手の receive_force
メソッドを呼び出します.その際に引数として力ベクトルを渡します.力を及ぼすべき相手としては PointMass
オブジェクトが想定されており,属性
pos
, vel
, radius
, mass
, viscous_damping
, restitution
を読み取れるとします.
PointMass
は update
メソッドによって,自身に及ぼす力を計算するとともに,受け取った力による自身の運動を計算します.
このルールの下で各 Actor を実装することで,メインプログラム側では:
for a in actor_list:
a.update()
actor_list[:] = [a for a in actor_list if a.is_alive]
のように書くだけで,a
の具体的な型 (PointMass
か Spring
かなど) に関わらずすべての Actor が協調して動きます.また:
for a in actor_list:
a.draw(screen)
のように書くことで a
の具体的な型に関わらず必要な描画が行われます.その「update」によってバネ復元力を発生するのか反発力を発生するのか運動計算をするのかを決めるは各 Actor の責務で,「draw」によって何が描かれるかも各 Actor の責務です.この性質をポリモーフィズムと呼ぶのでした.
Python だから PointMass
でも Spring
でも構わずに
actor_list
に入れられますが,Java や C++ のようにリストの型を指定しなくてはいけない言語ではそうはいかないと思います.
6章で Particle
と ConfinedParticle
を使い分けたときは,継承関係があったから基底クラスである Particle
型のリストを使えばよかったのですが,今回は PointMass
や Spring
などには継承関係がありません.こういうときはどうすればよいですか?
PointMass
や Spring
などの共通の親となるような基底クラスを作ればよいです.この章で Actor と呼んでいるものは,実は概念的にはそのような共通の基底クラスを表しています.
Actor 自身のインスタンスを生成することはありません.派生クラスを作ることのみを目的として用意するものです.このようなクラスを抽象基底クラス (abstract base class, ABC) と呼びます.
質問者も言っているように Python では抽象基底クラスを定義する必要はありません.しかし,定義しておくとプログラムの意味が明確になり,型チェックなどもかけられるようになります.興味がある人は abc モジュールについて調べてみてください.
なお,ここで:
actor_list[:] = [a for a in actor_list if a.is_alive]
の [:]
は必須です.ここを actor_list = ...
としてしまうと,
Boundary
型や CollisionResolver
型のオブジェクトに渡した actor_list
とは別のものを指してしまうからです.
7.4.2. 質点¶
では spring_mass.py
の内容をすべて削除するところから始めます.その空のファイルに以下を記載します.
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
以降で pygame.math.Vector2
を多用するため,4行目で PgVector
という別名で使えるようにしています (ちょっと不自然ですが,実はこれは次の章への布石です).
World
は説明不要だと思います.particle.py
の World
とは微妙に違うことにだけ注意してください.
CircleDrawer
と LineDrawer
は既に説明した通りです.
さらに,質点クラスを定義する以下のコードを追記します.
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
は,重力と粘性抵抗力を自身に及ぼす役割と,受け取った力に応じて位置と速度を更新する役割の 2 つを持っています.前者は update
メソッドの中から generate_force
メソッドを呼び出すことで行います (73-76行目).
ただし \(m\), \(d\), \(\boldsymbol{v}\) はそれぞれ質量,粘性抵抗係数,速度です.計算した力ベクトルを,自身の receive_force
メソッドを呼び出して渡します.
後者の役割のために,total_force
という合力を蓄積しておく属性を用意し,最初はゼロベクトルにしておきます (60行目).
receive_force
メソッドで受け取った力ベクトルはここに足しこんでいきます (71行目).
update
メソッドでは,move
メソッドを呼ぶことでその時点までの合力から運動を計算します (78-80行目).加速度が力から計算されている以外は
particle.py
と同じです.
運動が確定したら,属性 total_force
をゼロにリセットしておきます
(65行目).
このように receive_force
を介して力計算と運動計算を分けるのは回りくどく見えるかもしれません.しかしこうすることにより,バネや壁など他の要素から受け取る力と統一的に扱えるようになります.
残るメソッドは draw
です.初期化時に用意しておいた描画ストラテジ
self.drawer
を呼び出します.
self.drawer
を初期化しているのは 51 行目ですが,ここの
drawer or CircleDrawer("blue", 1)
という書き方は見慣れない人が多いと思います.これは, drawer
が None
でなければその値にし,
None
ならば or
の後ろの値にするという意味です.このように「もし値が未定ならばこの値にする」というときによく使う書き方です.
drawer or CircleDrawer("blue", 1)
の部分がなぜそのような動作になるのかわかりません.まず or
は真偽値に対して「または」の演算をすることを思い出してください.「x または y」は,x と y のどちらか一方が真であれば,もう一方の値をチェックすることなく真になります.
このことを利用して,Python を含む多くのプログラミング言語では,
or
演算の対象を左から順に評価していき,真であるものが見つかった時点でその値を返し,それより後のものは無視します.
and
も同様で,左から順に見ていって偽が現れた時点でその値を返します.(C 言語の ||
や &&
も同じです)
なぜそのような仕様になっているかというと,以下のような書き方をするときに便利だからです:
if k >= 0 and k < len(array) and array[k] == some_value:
do_something()
この例では,整数 k
がシーケンス array
のインデックスの範囲内であることを確認してから array[k]
の値を調べています.もし and
演算子が左から見ていくと保証されていなかったとすると,インデックス範囲をチェックする前に array[k]
を評価してエラーになってしまうかもしれません.
元の例では,Python の None
は真偽値として評価されたときに偽になることに注意します.もし drawer
が None
なら後ろの
CircleDrawer("blue", 1)
に評価が移りこれが結果になります.もし
drawer
に描画オブジェクトが渡されている場合は偽になりませんので,後ろの CircleDrawer("blue", 1)
は無視されて,
drawer
の値が結果になります.
このような書き方を自分でも使ってみようと思う人は,Python で偽と評価されるのは None
だけではないことに注意してください.例えば 0
,空文字列 ""
,空リスト []
なども真偽値として評価すると偽になりますので, or
演算にかけてもそこで止まらず次に進んでしまいます.
ここまでの動作を確認するため,
interactive_main.py
から CollisionResolver
型オブジェクトと
Boundary
型オブジェクトを生成している部分をいったんコメントアウトしておきます.
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
これで,左クリックで質点を生成することだけはできるようになったはずです.重力加速度や初速度はお好みで調整してください.
右クリック (固定質点の生成) やスペースを押しながらのクリック (バネ生成) をするとエラーになるので注意してください.まだ作ってないからです.
7.4.3. 固定質点¶
次に,PointMass
を継承して FixedPointMass
を作ります.以下のコードを
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
move
メソッドを上書きして動かなくする (90-91行目) だけで済みそうに思うのですが,実は速度と質量も上書きしておかないと衝突計算が正しくできません.
PointMass
の __init__
メソッドの内容をコピーペーストするのは無駄が多いので,PointMass
の __init__
メソッドを呼び出した上で,
self.vel
にゼロベクトルを,
self.mass
に非常に大きな値を代入し直すことにします.
super().__init__
は継承元 (スーパークラス) の __init__
メソッドを表す Python の記法です.
7.4.4. バネ¶
続いてバネを定義します.
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)
Spring
クラスの update
メソッドでは,まず両端の質点の存在を確認して,どちらかが存在しなければ自身の is_alive
を False
にします
(117-119行目).両方存在するなら generate_force
メソッドで力を及ぼします.
バネ復元力の計算自体は難しいことはないです.質点1が受け取る力は
です.ただし \(\boldsymbol{n}_\text{1,2}\) は質点1から2に向かう単位ベクトルで,\(k\), \(L\), \(L_\text{natural}\) はそれぞれバネ定数,現在の長さ,自然長です.
95-100行目は,基本的にはこの式の通りですが,注意が必要なのは単位ベクトルを求めるために距離で割っているところで (99行目),距離が 0 だとゼロ除算が発生します.質点が有限の大きさを持っていて (だからそれは質点じゃないだろというツッコミはスルーします) 衝突判定が有効なら本来この状況は生じませんが,シミュレーション上では生じ得ます.ここでは単純に,その場合は力を発生しないことにします.
interactive_main.py
を実行して,Space キーを押しながら質点を配置してバネの動作を確認してください.重力加速度や初速度をゼロにしておく方がわかりやすいかも知れません.
test_spring_mass.py
ver 14.0).test_spring_mass.py
を作成して以下を記載し,Explorer サイドバーから Run All Tests を実行してください.このフォルダではまだ pytest が有効化されていないので,有効化を促すポップアップが出るはずです.particle
フォルダと同じく有効化してください.
関数 generate_restoring_force
をテストする関数を作ります.バネ復元力は「正解の値」を用意するのが簡単なので,引数 f1_expected
として与えてやることにしています.また,None
が返ることが期待される場合 (復元力の計算しない場合) は引数 expects_none
を
True
にします.これらの引数は @pytest.mark.parametrize
で与えます.
x 方向と y 方向それぞれに伸びているバネしかテストしていません.本当は斜め方向もテストするべきですが,被テスト関数側ではベクトル表記で計算しているので,独立な 2 方向が正しく動けばまあ大丈夫かなということで省略しています.気になる人は,45度や60度などのわかりやすい方向のテストを加えてみるとよいと思います.
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)
意外に聞こえるかもしれませんが,まずテストを書き,テストに通るようにプログラム本体を作っていくというやり方が有効にはたらく場合があります.テスト駆動開発 (test driven development, TDD) と呼ばれます.
バネ復元力の例で考えてみましょう.一つ前の Q&A の単体テストの説明を読んでおいてください.
さて,自然長が 0 の場合であれば,特に難しいことは考えずに作れると思います.だからまずその場合に限定したテストのみを用意します.
x 方向に伸びる場合だけを考えましょう.期待される答え
f1_expected
は暗算で出せますね:
@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)
そして,このテストを通すような compute_restoring_force
を作るのも簡単です:
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
次に自然長が 0 でない場合のテストを追加していきます:
@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),
バネ長 = 10 よりも自然長 nl
が短い場合を用意しました.これらの場合に発生する力の符号がどちらになるかは,簡単にわかるので
f1_expected
として指定しておきます.
先ほど作った compute_restoring_force
では新しい 2 ケースのテストは通らないはずです.通るように直すのは難しくありません:
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
同様に,バネ長より自然長が長い場合のテストを追加します:
@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),
テストを走らせてみると,これら新しいテストケースも通るはずです.さっきのコードのままでよさそうです.
y 方向のテストも加えてみます.ついでにちょっとずれた位置なんかも試してみましょう:
@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),
特に問題なく通るようです.
ここまで来ると,基本的なロジックは大丈夫そうだなと自信が持てるようになります.気になる人は斜め方向のテストを加えてもよいでしょう.
あとは,例外的な条件が存在しないかを確認します.作成した
compute_restoring_force
を読み直すと,distance
で割っているところがあるので,ゼロ除算のおそれがあります.テストケースとして追加しましょう:
@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
これらのテストケースは通りません.距離が 0 のときは None
を返すように修正します:
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
これで ver 14.0 のコードにたどり着きました.
コツは,まず最初は例外的なケースや細かい検討を要するケースを無視して,問題のコアな部分を素直に実装することです.その後,細かいケースのテストを追加しながら精密化していくようにします.
7.4.5. 破断するバネ¶
Spring
を拡張する例として,あるしきい値を超える力がかかると破断するバネ FragileSpring
を作ってみます.
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
FragileSpring
の __init__
メソッドでは,スーパークラスの
__init__
を呼んだ上で,break_threshold
属性を追加しています
(137-139行目).
generate_force
メソッドでは,スーパークラスの generate_force
を呼び出して…としたいところなのですが,力の計算結果を得る手段がないので,仕方なく同じ処理内容をコピーペーストして作ります.その上で破断処理
(147-148行目) を追加します.
interactive_main.py
で,ファクトリの create_spring
メソッドが FragileSpring
を生成するようにしておきましょう.
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)
このくらいの定数設定だと,画面の端から端までの長さのバネを作ると消えてなくなるはずです.
7.4.6. 衝突解決¶
次は CollisionResolver
です.難しそうに感じますが,基本は他の Actor と同じで,力を及ぼすべき相手の receive_force メソッドを呼んで力ベクトルを渡していくのが責務です.
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)
175 行目で描画ストラテジ self.drawer
を初期化しています.引数 drawer
が None
でなければその値を,
None
なら「何も描画しない関数」をラムダ式で指定しています.
update
メソッドは generate_force
メソッドを呼び出すだけです
(184行目).generate_force
メソッドは,まず力を及ぼすべき相手をリストアップし (191行目),それらの相互間の衝突判定をし,衝突しているなら力を発生します.
191行目では,actor_list
属性に入っている Actor のうち対象となる条件を満たすものだけを変数 plist
に取り出します.条件は
target_condition
関数で判定しますが,デフォルト (178-179行目) では
PointMass
クラスのインスタンスかどうか (151-152行目) で判断します.条件は Strategy パターンで与えているので必要に応じて入れ替えられます.
ここで使っている isinstance(x, A)
は Python が用意する関数で,オブジェクト x
がクラス A
のインスタンスなら True
になります.ただし,A
自体ではなく A
のサブクラスのインスタンスでも True
になるのが特徴です.なので,FixedPointMass
も衝突計算の対象になります.
plist
に対象オブジェクト \(n\) 個がリストアップできたら,そのうち異なる 2 個を取り出す組合せを列挙します.一方を plist
の \(i\) 番目,他方を \(j\) 番目とすると,
\(i\) として \(0\) から \(n-1\) までを列挙しつつ,
\(j\) として \(i + 1\) から \(n-1\) まで列挙することで重複の無い全組合せが得られます (193-194行目).
衝突計算は compute_impact_force_between_points
関数で行います
(155-168行目).まず,中心間の距離が半径の和より大きければ衝突はしていませんので除外します (156-157行目).衝突している場合は,中心間を結ぶベクトルに沿った方向に弾性衝突による撃力が発生するとします.
さて,ここからしばらく高校物理の時間です.質点1から質点2に向かう方向の成分だけの 1 次元運動を考えます.質点1と質点2の質量と速度をそれぞれ \(m_1\), \(m_2\), \(v_1\), \(v_2\) とすると,運動量保存則は
と表され,弾性反発係数を \(e\) とすると
を満たす必要があります.そしてこれらを満たす速度変化を引き起こすために質点1と2に加えられるべき力は以下で与えられます.
これらを解くと以下が得られます.
というわけで,この式に従って,質点1から2に向かう方向の力を発生することになります.この方向の単位ベクトルを求め (160 行目),その方向の速度成分を pygame.math.Vector2
の dot
メソッドによる内積計算で求めます (161-162行目).あとはこの単位方向ベクトルに力の大きさを乗じれば OK
です (167行目).
ただし,いくつか考慮すべき細かい点があります.
- まず,質点1と2の位置が完全に一致していたら力の方向が定まりません.バネのときと同様にこれも本来はあり得ない状況ですが,計算上は起こり得ます.やはり単に無視することにします (158-159行目).
- 次に,
particle.py
で壁との衝突を考えたときにも発生した,衝突判定が複数回生じてしまう現象への対策も必要です.particle.py
と同じく速度条件を考えます.つまり両者が離れていく相対速度になっているときは除外します (163-164行目). - 最後に,反発係数をどのように与えるかです.質点1と質点2のそれぞれに属性
restitution
が設けられていますが,本来,反発係数は衝突する二体の組ごとに与えられるもので,質点単体に定義できるものではありません.しかし,衝突の可能性のある二体ごとに利用者が定義するのは面倒なので,多くの物理シミュレータでは,質点単体に反発係数 (っぽい何か) を定義し,そこから何らかのルールで計算に用いる反発係数を求めるという便宜的方法を採用しています.ここでは両者のrestitution
属性の積を使うことにします.物理的根拠はありません.
ここまでできたら,interactive_main.py
に CollisionResolver
を復帰させて動作を確認してください.
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"))
None
を返すのはなぜですか? PgVector((0, 0))
を返すようにすればもっとすっきりしたコードになりますよね?バネの場合はあまり問題にならないのですが,衝突の場合はチェックする組み合わせの数が多くなります (質点の数の2乗に比例します).そして通常ほとんどの場合に衝突は起きておらず,本来何もする必要がありません.PgVector((0, 0))
を返しても計算結果は同じですが,
receive_force
メソッドの呼び出しが無駄に頻繁に発生して遅くなります.
test_spring_mass.py
ver 16.0).バネ復元力より数式が複雑なので,個々の「正解」を外から与えてやるのは面倒です.かといって,被テスト関数
compute_impact_force_between_points
の中で計算している式をそのままテスト関数に書いても,テストの意義が弱いです.
ここではちょっと発想を変えて,被テスト関数が返す力によって変化する前後の速度を計算し,運動量が保存されていることと,弾性衝突の関係が成り立っていることを確認することにします.
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
クラスです.
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)
先ほどの CollisionResolver
と似た構造をしていますが少し簡単です.
generate_force
メソッドでは衝突対象を列挙し (238行目),それら1つずつに対して衝突計算をします.
境界面の正規化した法線ベクトルを \(\boldsymbol{n}\),境界面に含まれる 1 点の座標を \(\boldsymbol{p}_0\) とします.法線ベクトルは壁面の内部に向かう方向に取っていたことを思い出しておいてください.
このとき,質点の中心座標 \(\boldsymbol{x}\) が壁内への侵入量の法線方向成分は内積 \(\boldsymbol{n} \cdot (\boldsymbol{x} - \boldsymbol{p}_0)\) により求められます.これに半径分のオフセットを加えたものが正であり,かつ速度の法線方向成分が正のとき,衝突が生じたと判定します (204-205行目).
その際の撃力は,改めて計算してもよいのですが,質点間の衝突のときに求めた式で,相手方の速度を \(v_2 = 0\),質量を無限大とした極限 \(m_2 \rightarrow \infty\) を考えるのが手っ取り早いです (206-209行目).
interactive_main.py
で,壁や床を作ったり外したりして遊びましょう.
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
test_spring_mass.py
ver 17.0).こちらも被テスト関数が返す力をそのままテストしても面白くないので,力によって変化する前後の速度を計算し,弾性衝突の関係が成り立っていることを確認することにします.
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. 床面への潜り込みの回避¶
最後に一つ問題が残っています.particle.py
でも問題になった,床面への沈み込みです.
particle.py
のときと同様に y 座標を強制的に床面にすることで回避したいと思いますが,今回は少し微妙な問題が生じます.これは PointMass
の責務でしょうか.Boundary
の責務でしょうか.
質点の位置計算に関することなので,本来的には PointMass
の責務です.しかし,PointMass
型オブジェクトは自身が床と衝突したことを知る方法がありません.バネからの力も,他の質点からの反発力も,床からの反発力も,区別せずに receive_force
メソッドで受け取っているだけだからです.
では Boundary
が質点の y 座標を変更するべきかというと,オブジェクトの属性を外から変更するのは避けるという方針に反するのであまりやりたくありません.
spring_mass モジュールをさらに拡充していこうと思うと,同じような問題は他の場合にも起こり得ます.例えば,他の質点と一定回数衝突したら壊れる
FragilePointMass
を作りたいとしましょう.しかし他の質点との衝突によって力を受け取ったことを,それ以外の力の発生と区別する手段がありません.
要するに,generate_force
する側から receive_force
する側に,力以外の情報を伝達する手段が必要そうです.ここでは以下のような仕組みを作ります.
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})
- 61行目, 67行目
- 力を受け取る側の属性として,合力とは別にメッセージを保存するリスト
message_list
を用意します.このリストは最初は空に初期化され (61行目),update
されるたびに空にリセットされます (67行目). - 75-76行目
receive_force
メソッドとは別にreceive_message
メソッドを用意します.- 258-259行目
- メッセージを送る方を先に説明します.
Boundary
のgenerate_force
メソッドは,256行目で力を発生させた後,もし自身が床面オブジェクトだったら相手先にメッセージを送信します.メッセージは Python の辞書として作り,"type"
キーに対する値を"floor_hit"
,"y"
キーに対する値を y 座標とします. - 87-91行目
PointMass
に戻ります.メッセージの処理は,pygame のイベントの処理に倣って作ります.message_list
からメッセージを一つずつ取り出し,type が"floor_hit"
で,かつ自身が下向きの速度成分を持っていたら,メッセージの"y"
キーの値を使って自身の y 座標を更新します.
実装するのは簡単ではありませんが,基本的な考え方だけ説明しておきます.
質点は壁の内部に入り込まないという拘束条件を満たしているべきです.この拘束を満たさない場合のペナルティとなる値 (例えば侵入量の二乗) を定義し,運動計算の都度,その値を最小化するように位置の調整を行うという方法が考えられます.
同様の拘束条件は,壁との間だけではなく質点どうしの間でも満たされるべきです.(現在のプログラムはこれを全く考慮していません.複数の FixedPointMass
をすり鉢状に並べて,そこに PointMass
を落としてみたりするとよくわかります).
そのほかにも,2 つの質点が剛体リンクでつながっているとか,伸び縮みしない紐で結ばれているといった状況も,質点の間の距離に関する拘束条件として表現できます.
このように複数の運動体の間に複数の拘束条件が存在する場合は,全体のペナルティが最小になるように質点系全体の位置を最適化することが必要になります.例えば非線形最小二乗問題として定式化されます.
Spring
や CollisionResolver
などと比べると,なんだか
PointMass
だけ構造が妙に複雑なように思うのですが?PointMass
は抱えている責務が多すぎます.本来は複数のクラスに分割すべきです.PointMass
は力発生と運動計算という 2 つの責務を同時に担っています.物理学的に見ても質量には重力質量と慣性質量という 2 つの側面があり,本来は別物なのに一致しているものと理解されます.
少なくともこれら 2 つの責務は別のクラス (例えば GravityGenerator
と MovingPoint
など) に分け,PointMass
クラスはそれぞれのオブジェクトへの参照を保持して (オブジェクト合成),処理を振り分けるのが望ましい姿だと思います.
別の方針として,PointMass
をこれら 2 つのクラスを継承することにより定義することも考えられます.しかし,複数のクラスを継承する多重継承(multiple inheritance) は一般に設計が難しく,下手に手を出さない方がよいとされています.Java のように多重継承を禁止している言語もあります.
7.5. 演習¶
例題7-1
独自の LineDrawer を定義して,Spring
の表示が点線になるようにしてください.例えば 2 つの端点を結ぶ線分を適当な個数の区間に分割し,各区間につき 1 個の円を表示するなどの方法を考えるとよいです.
例題7-2
重力の無い世界で,バネ定数 \(k\) のバネの一端を固定し,もう一端につないだ質量 \(m\) の質点に,正弦波の強制振動力を加えることを考えます.強制振動の角周波数が共振角周波数 \(\sqrt{k / m}\) に近い場合と遠い場合の挙動をシミュレーションで確認してください.
(ヒント: 角周波数 \(\omega\) の単振動は,フレーム番号を \(n\), フレームあたりの時間ステップを \(\Delta t\) とすると \(\cos(\omega n \Delta t)\) で表せます.振幅は適宜設定してください)
例題7-3
CollisionResolver
を継承し,累計衝突回数を画面の左上に数字で表示する機能を持つ衝突解決クラスを作成してください.
(ヒント: 例えば,累計衝突回数を保持する属性を追加し,
generate_force
を上書きしてそれをカウントアップするとともに,
draw
を上書きしてそれを表示するといった方法が考えられます)
例題7-4
Spring
クラスの実装を参照しながら,2 つの質点の間に距離の二乗に反比例する引力や斥力を発生させる Actor のクラスを作成してください.