0_getting_started

Author: liminchen

0_getting_started/readme.md

@@ -0,0 +1,22 @@
+# Getting Started with Solid Simulation
+
+Free fall, spring, and mass-spring solid simulation with symplectic Euler time integration.
+
+## Dependencies
+```
+pip install numpy scipy pygame
+```
+
+## Run
+Free falling particle:
+```
+python simulator0.py
+```
+Hanging particle:
+```
+python simulator1.py
+```
+Hanging square:
+```
+python simulator2.py
+```

0_getting_started/simulator0.py

@@ -0,0 +1,52 @@
+# Free Fall Simulation
+
+import math
+import numpy as np  # for vector data structure and computations
+import pygame       # for visualization
+
+# simulation setup
+m = 1000                    # mass of particle
+x = np.array([0.0, 0.0])    # position of particle
+v = np.array([0.1, 0.1])    # velocity of particle
+g = np.array([0, -9.81])    # gravitational acceleration
+h = 0.1                     # time step size in seconds
+
+# visualization/rendering setup
+pygame.init()
+render_FPS = 100                    # number of frames to render per second
+resolution = np.array([900, 900])   # visualization window size in pixels
+offset = resolution / 2             # offset between window coordinates and simulated coordinates
+scale = 200                         # scale between window coordinates and simulated coordinates
+def screen_projection(x):           # convert simulated coordinates to window coordinates
+    return [offset[0] + scale * x[0], resolution[1] - (offset[1] + scale * x[1])]
+screen = pygame.display.set_mode(resolution)    # initialize visualizer
+
+time_step = 0   # the number of the current time step
+running = True  # flag indicating whether the simulation is still running
+while running:
+    # run until the user asks to quit
+    for event in pygame.event.get():
+        if event.type == pygame.QUIT:
+            running = False
+    
+    # update the frame to display according to render_FPS
+    if time_step % int(math.ceil((1.0 / render_FPS) / h)) == 0:
+        # fill the background with white color, display simulation time at the top,
+        # render a floor at y=-1, and render the particle as a circle:
+        screen.fill((255, 255, 255))    
+        pygame.display.set_caption('Current time: ' + f'{time_step * h: .2f}s')
+        pygame.draw.aaline(screen, (0, 0, 255), screen_projection([-2, -1]), screen_projection([2, -1]))
+        pygame.draw.circle(screen, (0, 0, 255), screen_projection(x), 0.1 * scale)
+        pygame.display.flip()   # flip the display
+        pygame.time.wait(int(1000.0 / render_FPS))  # wait to render the next frame
+
+    # step forward the simulation by updating particle velocity and position
+    v += h * g / m
+    x += h * v
+
+    # pause the simulation when the particle touches on the ground 
+    if x[1] <= -1:
+        input()
+        break
+
+    time_step += 1  # update time step counter

0_getting_started/simulator1.py

@@ -0,0 +1,58 @@
+# Spring Simulation
+
+import math
+import numpy as np  # for vector data structure and computations
+import pygame       # for visualization
+
+# simulation setup
+m = 1000                    # mass of particle
+x = np.array([0.1, 0.0])    # position of particle
+v = np.array([0.0, 0.0])    # velocity of particle
+g = np.array([0, -9.81])    # gravitational acceleration
+spring_rest_len = 0.3       # rest length of the spring ###
+spring_stiffness = 1000     # stiffness of the spring ###
+h = 0.1                     # time step size in seconds
+
+# visualization/rendering setup
+pygame.init()
+render_FPS = 100                    # number of frames to render per second
+resolution = np.array([900, 900])   # visualization window size in pixels
+offset = resolution / 2             # offset between window coordinates and simulated coordinates
+scale = 200                         # scale between window coordinates and simulated coordinates
+def screen_projection(x):           # convert simulated coordinates to window coordinates
+    return [offset[0] + scale * x[0], resolution[1] - (offset[1] + scale * x[1])]
+screen = pygame.display.set_mode(resolution)    # initialize visualizer
+
+time_step = 0   # the number of the current time step
+running = True  # flag indicating whether the simulation is still running
+while running:
+    # run until the user asks to quit
+    for event in pygame.event.get():
+        if event.type == pygame.QUIT:
+            running = False
+    
+    # update the frame to display according to render_FPS
+    if time_step % int(math.ceil((1.0 / render_FPS) / h)) == 0:
+        # fill the background with white color, display simulation time at the top,
+        # render a floor at y=-1, draw the spring segment, and render the particle as a circle:
+        screen.fill((255, 255, 255))    
+        pygame.display.set_caption('Current time: ' + f'{time_step * h: .2f}s')
+        pygame.draw.aaline(screen, (0, 0, 255), screen_projection([-2, -1]), screen_projection([2, -1]))
+        pygame.draw.aaline(screen, (0, 0, 255), screen_projection([0, 0]), screen_projection(x)) ###
+        pygame.draw.circle(screen, (0, 0, 255), screen_projection(x), 0.1 * scale)
+        pygame.display.flip()   # flip the display
+        pygame.time.wait(int(1000.0 / render_FPS))  # wait to render the next frame
+
+    # step forward the simulation by updating particle velocity and position
+    spring_cur_len = math.sqrt(x[0] * x[0] + x[1] * x[1]) ###
+    spring_displacement = spring_cur_len - spring_rest_len ###
+    spring_force = -spring_stiffness * spring_displacement * (x / spring_cur_len) ###
+    v += h * (g + spring_force) / m
+    x += h * v
+
+    # pause the simulation when the particle touches on the ground 
+    if x[1] <= -1:
+        input()
+        break
+
+    time_step += 1  # update time step counter

0_getting_started/simulator2.py

@@ -0,0 +1,69 @@
+# Mass-Spring Solid Simulation
+
+import math
+import numpy as np  # for vector data structure and computations
+import pygame       # for visualization
+import square_mesh  # for generating a square mesh
+
+# simulation setup
+m = 1000                    # mass of each particle
+side_length = 1             # side length of the square
+n_seg = 4                   # number of springs per side of the square
+[x, e] = square_mesh.generate(side_length, n_seg)   # array of particle positions and springs   ###
+v = np.array([[0.0, 0.0]] * len(x))     # velocity array of particles ###
+g = np.array([0, -9.81])    # gravitational acceleration
+spring_rest_len = []        # rest length array of the springs ###
+for i in range(0, len(e)):  # calculate the rest length of each spring
+    spring_vec = x[e[i][0]] - x[e[i][1]]    # the vector connecting two ends of spring i
+    spring_rest_len.append(math.sqrt(spring_vec[0] * spring_vec[0] + spring_vec[1] * spring_vec[1]))
+spring_stiffness = 1000     # stiffness of the spring
+h = 0.1                     # time step size in seconds
+
+# visualization/rendering setup
+pygame.init()
+render_FPS = 100                    # number of frames to render per second
+resolution = np.array([900, 900])   # visualization window size in pixels
+offset = resolution / 2             # offset between window coordinates and simulated coordinates
+scale = 200                         # scale between window coordinates and simulated coordinates
+def screen_projection(x):           # convert simulated coordinates to window coordinates
+    return [offset[0] + scale * x[0], resolution[1] - (offset[1] + scale * x[1])]
+screen = pygame.display.set_mode(resolution)    # initialize visualizer
+
+time_step = 0   # the number of the current time step
+running = True  # flag indicating whether the simulation is still running
+while running:
+    # run until the user asks to quit
+    for event in pygame.event.get():
+        if event.type == pygame.QUIT:
+            running = False
+    
+    # update the frame to display according to render_FPS
+    if time_step % int(math.ceil((1.0 / render_FPS) / h)) == 0:
+        # fill the background with white color, display simulation time at the top,
+        # draw each spring segment, and render each particle as a circle:
+        screen.fill((255, 255, 255))    
+        pygame.display.set_caption('Current time: ' + f'{time_step * h: .2f}s')
+        for i in range(0, len(e)):  ###
+            pygame.draw.aaline(screen, (0, 0, 255), screen_projection(x[e[i][0]]), screen_projection(x[e[i][1]]))
+        for i in range(0, len(x)):  ###
+            pygame.draw.circle(screen, (0, 0, 255), screen_projection(x[i]), 0.02 * scale)
+        pygame.display.flip()   # flip the display
+        pygame.time.wait(int(1000.0 / render_FPS))  # wait to render the next frame
+
+    # step forward the simulation by updating particle velocity and position ###
+    for i in range(0, len(e)):
+        # calculate elasticity force of spring i:
+        spring_vec = x[e[i][0]] - x[e[i][1]]
+        spring_cur_len = math.sqrt(spring_vec[0] * spring_vec[0] + spring_vec[1] * spring_vec[1])
+        spring_displacement = spring_cur_len - spring_rest_len[i]
+        spring_force = -spring_stiffness * spring_displacement * (spring_vec / spring_cur_len)
+        # update the velocity of the two ends of spring i
+        v[e[i][0]] += h * (g + spring_force) / m
+        v[e[i][1]] += h * (g - spring_force) / m
+    # fix the top left and top right corner by setting velocity to 0:
+    v[n_seg] = v[(n_seg + 1) * (n_seg + 1) - 1] = np.array([0, 0])
+    # update the position of each particle:
+    for i in range(0, len(x)):
+        x[i] += h * v[i]
+
+    time_step += 1  # update time step counter

0_getting_started/square_mesh.py

@@ -0,0 +1,27 @@
+import numpy as np
+
+def generate(side_length, n_seg):
+    # sample nodes uniformly on a square
+    x = np.array([[0.0, 0.0]] * ((n_seg + 1) ** 2))
+    step = side_length / n_seg
+    for i in range(0, n_seg + 1):
+        for j in range(0, n_seg + 1):
+            x[i * (n_seg + 1) + j] = [-side_length / 2 + i * step, -side_length / 2 + j * step]
+    
+    # connect the nodes with edges
+    e = []
+    # horizontal edges
+    for i in range(0, n_seg):
+        for j in range(0, n_seg + 1):
+            e.append([i * (n_seg + 1) + j, (i + 1) * (n_seg + 1) + j])
+    # vertical edges
+    for i in range(0, n_seg + 1):
+        for j in range(0, n_seg):
+            e.append([i * (n_seg + 1) + j, i * (n_seg + 1) + j + 1])
+    # diagonals
+    for i in range(0, n_seg):
+        for j in range(0, n_seg):
+            e.append([i * (n_seg + 1) + j, (i + 1) * (n_seg + 1) + j + 1])
+            e.append([(i + 1) * (n_seg + 1) + j, i * (n_seg + 1) + j + 1])
+
+    return [x, e]