steffen/VoxelEngine
1
0
Fork 0
Code Issues Pull requests Releases Wiki Activity

fluid simu kinda

Browse source
This commit is contained in:
zomseffen 2020-11-30 19:15:08 +01:00
parent 41761b7b56
commit 6577ac6558
1 changed files with 158 additions and 151 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

301
Client/Client.py
View file

@ -21,9 +21,10 @@ import time
from scipy.signal import convolve from scipy.signal import convolve
MAX_DISTANCE = 200.0 MAX_DISTANCE = 200.0
FRICTION_COEFFICENT = 0.9 FRICTION_COEFFICENT = 1
EPSILON = 0.00001 EPSILON = 0.00001
def value_to_color(v, min_value, max_value): def value_to_color(v, min_value, max_value):
r = g = b = 0.0 r = g = b = 0.0
scope = max_value - min_value scope = max_value - min_value
@ -38,6 +39,7 @@ def value_to_color(v, min_value, max_value):
b /= l b /= l
return r, g, b return r, g, b
class Client: class Client:
def __init__(self, test=False, pos=[0, 0, 0]): def __init__(self, test=False, pos=[0, 0, 0]):
with open('./config.json', 'r') as f: with open('./config.json', 'r') as f:
@ -90,7 +92,8 @@ class Client:
glAttachShader(self.normal_program[key], self.fragment_shader_id) glAttachShader(self.normal_program[key], self.fragment_shader_id)
glLinkProgram(self.normal_program[key]) glLinkProgram(self.normal_program[key])
self.depth_program[self.normal_program[key]] = Spotlight.getDepthProgram(self.vertex_shader_id, key.GeometryShaderId) self.depth_program[self.normal_program[key]] = Spotlight.getDepthProgram(self.vertex_shader_id,
key.GeometryShaderId)
self.world_provider = WorldProvider(self.normal_program) self.world_provider = WorldProvider(self.normal_program)
for x_pos in range(0, 100): for x_pos in range(0, 100):
@ -113,11 +116,44 @@ class Client:
self.time = time.time() self.time = time.time()
self.heat_map = np.zeros((100, 100, 1)) self.field = (100, 100, 1)
self.e_a = np.array([
[0, 0, 0],
[1, 0, 0],
[1, 1, 0],
[0, 1, 0],
[-1, 1, 0],
[-1, 0, 0],
[-1, -1, 0],
[0, -1, 0],
[1, -1, 0],
])
self.v_map_x = np.zeros((100, 100, 1)) self.relaxation_time = 0.55 # 0.55
self.v_map_y = np.zeros((100, 100, 1)) self.w_a = [
self.v_map_z = np.zeros((100, 100, 1)) 4.0 / 9.0,
1.0 / 9.0,
1.0 / 36.0,
1.0 / 9.0,
1.0 / 36.0,
1.0 / 9.0,
1.0 / 36.0,
1.0 / 9.0,
1.0 / 36.0
]
self.n_a = np.zeros((len(self.e_a),) + self.field)
self.n_a_eq = np.zeros(self.n_a.shape)
self.n = np.zeros(self.field)
self.n[:, :, :] += 1.0
# self.n /= np.sum(self.n)
self.n_a[0] = np.array(self.n)
self.u = np.zeros(self.field + (self.e_a.shape[1],))
self.compressible = True
self.max_n = self.w_a[0]
self.test_pixel = [40, 50, 0]
if not test: if not test:
glutMainLoop() glutMainLoop()
@ -125,8 +161,6 @@ class Client:
self.display() self.display()
self.resize(100, 100) self.resize(100, 100)
def display(self): def display(self):
glClearColor(0, 0, 0, 0) glClearColor(0, 0, 0, 0)
glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT) glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT)
@ -172,171 +206,143 @@ class Client:
glUniform1f(widthid, self.width) glUniform1f(widthid, self.width)
glUniform1f(heightid, self.height) glUniform1f(heightid, self.height)
world.render(translate(self.pos[0], self.pos[1], self.pos[2]) * lookAt(0, 0, 0, 0, 0, -self.pos[2], 0, 1, 0) * projMatrix, rotate(0, 0, 0)) world.render(translate(self.pos[0], self.pos[1], self.pos[2]) * lookAt(0, 0, 0, 0, 0, -self.pos[2], 0, 1,
0) * projMatrix, rotate(0, 0, 0))
glFlush() glFlush()
glutSwapBuffers() glutSwapBuffers()
max_value = np.max(self.heat_map)
# min_value = np.min(self.heat_map)
min_value = 0 min_value = 0
max_value_n = np.max(self.n)
# max_value = 1.0
vel = np.sqrt(np.square(self.v_map_x) + np.square(self.v_map_y) + np.square(self.v_map_z)) vel = np.sqrt(np.sum(np.square(self.u), axis=3))
max_value = np.max(vel) max_value_vel = np.max(vel)
min_value = np.min(vel) # max_value_vel = np.sqrt(3)
print('round')
print(max_value_n)
for x_pos in range(0, 100): for x_pos in range(0, 100):
for y_pos in range(0, 100): for y_pos in range(0, 100):
for z_pos in range(0, 1): for z_pos in range(0, 1):
# r, g, b = value_to_color(self.heat_map[x_pos, y_pos, z_pos], min_value, max_value) r, g, b = value_to_color(self.n[x_pos, y_pos, z_pos], min_value, max_value_n)
r, g, b = value_to_color(vel[x_pos, y_pos, z_pos], min_value, max_value) # r, g, b = value_to_color(vel[x_pos, y_pos, z_pos], min_value, max_value_vel)
self.world_provider.world.set_color(x_pos, y_pos, z_pos, r, g, b) self.world_provider.world.set_color(x_pos, y_pos, z_pos, r, g, b)
# friction self.world_provider.world.set_color(int(round(self.test_pixel[0])),
# self.heat_map += np.sqrt(np.square(self.v_map_x * (1.0 - FRICTION_COEFFICENT)) + int(round(self.test_pixel[1])),
# np.square(self.v_map_y * (1.0 - FRICTION_COEFFICENT)) + int(round(self.test_pixel[2])), 1.0, 1.0, 1.0)
# np.square(self.v_map_z * (1.0 - FRICTION_COEFFICENT)))
self.v_map_x *= FRICTION_COEFFICENT
self.v_map_y *= FRICTION_COEFFICENT
self.v_map_z *= FRICTION_COEFFICENT
# hot stuff rises / cool stuff sinks
rise = self.heat_map[:, :-1, :] > (self.heat_map[:, 1:, :] + EPSILON)
self.v_map_y[:, :-1, :] += 1.0 * rise
sink = self.heat_map[:, :-1, :] < (self.heat_map[:, 1:, :] - EPSILON)
self.v_map_y[:, 1:, :] -= 1.0 * sink
#flow
new_v_map_x = np.zeros(self.v_map_x.shape)
new_v_map_y = np.zeros(self.v_map_x.shape)
new_v_map_z = np.zeros(self.v_map_x.shape)
for x_pos in range(self.v_map_x.shape[0]):
for y_pos in range(self.v_map_x.shape[1]):
for z_pos in range(self.v_map_x.shape[2]):
target_x = min(self.v_map_x.shape[0] - 1,
max(0, int(round(x_pos + self.v_map_x[x_pos, y_pos, z_pos]))))
target_y = min(self.v_map_x.shape[1] - 1,
max(0, int(round(y_pos + self.v_map_y[x_pos, y_pos, z_pos]))))
target_z = min(self.v_map_x.shape[2] - 1,
max(0, int(round(z_pos + self.v_map_z[x_pos, y_pos, z_pos]))))
friction_dispersion = (1.0 -FRICTION_COEFFICENT) / 4 # self.u *= 0.95
# velocity dispersion
# x
new_v_map_x[target_x, target_y, target_z] += self.v_map_x[x_pos, y_pos, z_pos] * FRICTION_COEFFICENT
if target_y + 1 < self.v_map_x.shape[1] - 1:
new_v_map_y[target_x, target_y + 1, target_z] += self.v_map_x[x_pos, y_pos, z_pos] * friction_dispersion
else:
new_v_map_x[target_x, target_y, target_z] += self.v_map_x[x_pos, y_pos, z_pos] * friction_dispersion
if target_y - 1 > 0:
new_v_map_y[target_x, target_y - 1, target_z] -= self.v_map_x[x_pos, y_pos, z_pos] * friction_dispersion
else:
new_v_map_x[target_x, target_y, target_z] += self.v_map_x[x_pos, y_pos, z_pos] * friction_dispersion
if target_z + 1 < self.v_map_x.shape[2] - 1: old_n_sum = np.sum(self.n)
new_v_map_z[target_x, target_y, target_z + 1] += self.v_map_x[x_pos, y_pos, z_pos] * friction_dispersion for a in range(len(self.e_a)):
else: e_au = np.sum(self.e_a[a] * self.u, axis=3)
new_v_map_x[target_x, target_y, target_z] += self.v_map_x[x_pos, y_pos, z_pos] * friction_dispersion uu = np.sum(self.u * self.u, axis=3)
if target_z - 1 > 0: self.n_a_eq[a] = self.w_a[a] * self.n * (1.0 + 3.0 * e_au + 4.5 * np.square(e_au) - 1.5 * uu)
new_v_map_z[target_x, target_y, target_z - 1] -= self.v_map_x[x_pos, y_pos, z_pos] * friction_dispersion print(np.max(self.n_a_eq[0]))
else: if not self.compressible:
new_v_map_x[target_x, target_y, target_z] += self.v_map_x[x_pos, y_pos, z_pos] * friction_dispersion excess = (self.n_a_eq[0] > self.max_n) * (self.n_a_eq[0] - self.max_n)
# y dir_sum = np.sum(self.n_a_eq[1:], axis=0)
new_v_map_y[target_x, target_y, target_z] += self.v_map_y[x_pos, y_pos, z_pos] * FRICTION_COEFFICENT self.n_a_eq[1:] += excess * 1/8 #(self.n_a_eq[1:] / dir_sum)
if target_x + 1 < self.v_map_x.shape[0] - 1: self.n_a_eq[0] -= excess
new_v_map_x[target_x + 1, target_y, target_z] += self.v_map_y[x_pos, y_pos, z_pos] * friction_dispersion
else:
new_v_map_y[target_x, target_y, target_z] += self.v_map_y[x_pos, y_pos, z_pos] * friction_dispersion
if target_x - 1 > 0:
new_v_map_x[target_x - 1, target_y, target_z] -= self.v_map_y[x_pos, y_pos, z_pos] * friction_dispersion
else:
new_v_map_y[target_x, target_y, target_z] += self.v_map_y[x_pos, y_pos, z_pos] * friction_dispersion
if target_z + 1 < self.v_map_x.shape[2] - 1: n_a_t_1 = np.zeros((len(self.e_a),) + self.field)
new_v_map_z[target_x, target_y, target_z + 1] += self.v_map_y[x_pos, y_pos, z_pos] * friction_dispersion for a in range(len(self.e_a)):
temp = ((-1.0 / self.relaxation_time) * (self.n_a[a] - self.n_a_eq[a]) + self.n_a[a])
n_a_t_1[a, max(self.e_a[a][0], 0): min(self.field[0], self.field[0] + self.e_a[a][0]),
max(self.e_a[a][1], 0): min(self.field[1], self.field[1] + self.e_a[a][1]),
max(self.e_a[a][2], 0): min(self.field[2], self.field[2] + self.e_a[a][2])] += temp[
max(-self.e_a[a][0], 0): min(
self.field[0],
self.field[0] -
self.e_a[a][0]),
max(-self.e_a[a][1], 0): min(
self.field[1],
self.field[1] -
self.e_a[a][1]),
max(-self.e_a[a][2], 0): min(
self.field[2],
self.field[2] -
self.e_a[a][2])]
for index in range(len(self.e_a[a])):
if self.e_a[a][index] != 0:
e_a_clipped = -np.array(self.e_a[a])
# e_a_clipped[index] = 0
# e_a_clipped[index] = -self.e_a[a][index]
clipped_dir = np.zeros(3)
clipped_dir[index] = self.e_a[a][index]
e_a_index = np.where(np.all(self.e_a == e_a_clipped, axis=1))[0][0]
# e_a_index = 0
if index == 0:
if self.e_a[a][index] > 0:
slice_index = -1
else: else:
new_v_map_y[target_x, target_y, target_z] += self.v_map_y[x_pos, y_pos, z_pos] * friction_dispersion slice_index = 0
if target_z - 1 > 0: n_a_t_1[e_a_index, slice_index, :, :] += temp[slice_index, :, :]
new_v_map_z[target_x, target_y, target_z - 1] -= self.v_map_y[x_pos, y_pos, z_pos] * friction_dispersion temp[slice_index, :, :] *= 0
if index == 1:
if self.e_a[a][index] > 0:
slice_index = -1
else: else:
new_v_map_y[target_x, target_y, target_z] += self.v_map_y[x_pos, y_pos, z_pos] * friction_dispersion slice_index = 0
# z n_a_t_1[e_a_index, :, slice_index, :] += temp[:, slice_index, :]
new_v_map_z[target_x, target_y, target_z] += self.v_map_z[x_pos, y_pos, z_pos] * FRICTION_COEFFICENT temp[:, slice_index, :] *= 0
if target_x + 1 < self.v_map_x.shape[0] - 1:
new_v_map_x[target_x + 1, target_y, target_z] += self.v_map_z[x_pos, y_pos, z_pos] * friction_dispersion if index == 2:
if self.e_a[a][index] > 0:
slice_index = -1
else: else:
new_v_map_z[target_x, target_y, target_z] += self.v_map_z[x_pos, y_pos, z_pos] * friction_dispersion slice_index = 0
if target_x - 1 > 0: n_a_t_1[e_a_index, :, :, slice_index] += temp[:, :, slice_index]
new_v_map_x[target_x - 1, target_y, target_z] -= self.v_map_z[x_pos, y_pos, z_pos] * friction_dispersion temp[:, :, slice_index] *= 0
else:
new_v_map_z[target_x, target_y, target_z] += self.v_map_z[x_pos, y_pos, z_pos] * friction_dispersion
if target_y + 1 < self.v_map_x.shape[1] - 1: self.n_a = n_a_t_1
new_v_map_y[target_x, target_y + 1, target_z] += self.v_map_z[x_pos, y_pos, z_pos] * friction_dispersion if np.min(self.n_a) < 0:
else: test = 1
new_v_map_z[target_x, target_y, target_z] += self.v_map_z[x_pos, y_pos, z_pos] * friction_dispersion self.n_a = np.maximum(0, self.n_a)
if target_y - 1 > 0: self.n_a = old_n_sum * self.n_a / np.sum(self.n_a)
new_v_map_y[target_x, target_y - 1, target_z] -= self.v_map_z[x_pos, y_pos, z_pos] * friction_dispersion self.n = np.sum(self.n_a, axis=0, keepdims=False)
else: # self.n = np.sum(np.abs(n_a_t_1), axis=0, keepdims=False)
new_v_map_z[target_x, target_y, target_z] += self.v_map_z[x_pos, y_pos, z_pos] * friction_dispersion
# handle boundaries
filter_mat = np.array([[-1.0], [0], [1.0]]) / 2.0
new_v_map_y[0, :, :] += convolve(new_v_map_x[0, :, :], filter_mat, 'same')
new_v_map_x[0, :, :] = 0
new_v_map_y[new_v_map_x.shape[0] - 1, :, :] +=\
convolve(new_v_map_x[new_v_map_x.shape[0] - 1, :, :], filter_mat, 'same')
new_v_map_x[new_v_map_x.shape[0] - 1, :, :] = 0
filter_mat = np.array([[-1.0], [0], [1.0]]) / 2.0 # stabilise the number because of rounding errors
new_v_map_x[:, 0, :] += convolve(new_v_map_y[:, 0, :], filter_mat, 'same') self.n = old_n_sum * self.n / np.sum(self.n)
new_v_map_y[:, 0, :] = 0
new_v_map_x[:, new_v_map_x.shape[1] - 1, :] +=\
convolve(new_v_map_y[:, new_v_map_x.shape[1] - 1, :], filter_mat, 'same')
new_v_map_y[:, new_v_map_x.shape[1] - 1, :] = 0
# corners self.u *= 0
new_v_map_x[0, 0, 0] = new_v_map_y[0, 0, 0] = new_v_map_z[0, 0, 0] = 0 for a in range(len(self.e_a)):
new_v_map_x[-1, 0, 0] = new_v_map_y[-1, 0, 0] = new_v_map_z[-1, 0, 0] = 0 self.u[:, :, :, 0] += self.n_a[a] * self.e_a[a][0]
new_v_map_x[-1, -1, 0] = new_v_map_y[-1, -1, 0] = new_v_map_z[-1, -1, 0] = 0 self.u[:, :, :, 1] += self.n_a[a] * self.e_a[a][1]
new_v_map_x[-1, -1, -1] = new_v_map_y[-1, -1, -1] = new_v_map_z[-1, -1, -1] = 0 self.u[:, :, :, 2] += self.n_a[a] * self.e_a[a][2]
new_v_map_x[0, -1, -1] = new_v_map_y[0, -1, -1] = new_v_map_z[0, -1, -1] = 0 self.u[:, :, :, 0] /= self.n
new_v_map_x[0, -1, 0] = new_v_map_y[0, -1, 0] = new_v_map_z[0, -1, 0] = 0 self.u[:, :, :, 1] /= self.n
new_v_map_x[-1, -1, 0] = new_v_map_y[-1, -1, 0] = new_v_map_z[-1, -1, 0] = 0 self.u[:, :, :, 2] /= self.n
new_v_map_x[-1, 0, -1] = new_v_map_y[-1, 0, -1] = new_v_map_z[-1, 0, -1] = 0
self.v_map_x = new_v_map_x self.u[self.n == 0] = 0
self.v_map_y = new_v_map_y
self.v_map_z = new_v_map_z
filter_mat = (np.zeros((3, 3, 1)) + 1.0) / 9.0 length = np.sqrt(np.sum(np.square(self.u), axis=3, keepdims=True))
v_map = np.pad(self.v_map_x, 1, 'edge') # gravity
self.v_map_x = convolve(v_map, filter_mat, mode='same')[1:-1, 1:-1, 1:2] gravity_applies = self.n < self.w_a[0]
v_map = np.pad(self.v_map_y, 1, 'edge') gravity_applies[:, :-1, :] = gravity_applies[:, 1:, :]
self.v_map_y = convolve(v_map, filter_mat, mode='same')[1:-1, 1:-1, 1:2] gravity_applies[:, -1, :] = False
v_map = np.pad(self.v_map_z, 1, 'edge') self.u[gravity_applies, 1] -= 0.01
self.v_map_z = convolve(v_map, filter_mat, mode='same')[1:-1, 1:-1, 1:2]
# moving heat
heat_map = np.zeros(self.heat_map.shape)
for x_pos in range(self.v_map_x.shape[0]):
for y_pos in range(self.v_map_x.shape[1]):
for z_pos in range(self.v_map_x.shape[2]):
target_x = min(self.v_map_x.shape[0] - 1,
max(0, int(round(x_pos + self.v_map_x[x_pos, y_pos, z_pos]))))
target_y = min(self.v_map_x.shape[1] - 1,
max(0, int(round(y_pos + self.v_map_y[x_pos, y_pos, z_pos]))))
target_z = min(self.v_map_x.shape[2] - 1,
max(0, int(round(z_pos + self.v_map_z[x_pos, y_pos, z_pos]))))
heat_map[target_x, target_y, target_z] += self.heat_map[x_pos, y_pos, z_pos]
self.heat_map = heat_map new_lengths = np.sqrt(np.sum(np.square(self.u), axis=3, keepdims=True))
# dispersing heat self.u = self.u / new_lengths * length
heat_map = np.pad(self.heat_map, 1, 'edge') zero_length = (new_lengths == 0)
self.heat_map = convolve(heat_map, filter_mat, mode='same')[1:-1, 1:-1, 1:2] self.u[zero_length[:, :, :, 0], :] = 0
# heat source keeps source block on constant heat
self.heat_map[50-5:50+5, 0, 0] = 100.0
# roof gets cooled off to min temp
self.heat_map[:, 99, :] = np.maximum(self.heat_map[:, 99, :] * 0.8, 0.0)
print(1.0 / (time.time() - self.time)) u = self.u[int(self.test_pixel[0]), int(self.test_pixel[1]), int(self.test_pixel[2])]
self.test_pixel[0] = max(0, min(self.field[0] - 1, self.test_pixel[0] + u[0]))
self.test_pixel[1] = max(0, min(self.field[1] - 1, self.test_pixel[1] + u[1]))
self.test_pixel[2] = max(0, min(self.field[2] - 1, self.test_pixel[2] + u[2]))
# print(1.0 / (time.time() - self.time))
self.time = time.time() self.time = time.time()
glutPostRedisplay() glutPostRedisplay()
@ -382,5 +388,6 @@ class Client:
glutFullScreenToggle() glutFullScreenToggle()
# print(key) # print(key)
if __name__ == '__main__': if __name__ == '__main__':
client = Client(pos=[-50, -50, -200]) client = Client(pos=[-50, -50, -200])