From 6577ac6558873101a8c2c3df5680e25bf1cf505e Mon Sep 17 00:00:00 2001 From: zomseffen Date: Mon, 30 Nov 2020 19:15:08 +0100 Subject: [PATCH] fluid simu kinda --- Client/Client.py | 309 ++++++++++++++++++++++++----------------------- 1 file changed, 158 insertions(+), 151 deletions(-) diff --git a/Client/Client.py b/Client/Client.py index dfe456e..e36563a 100644 --- a/Client/Client.py +++ b/Client/Client.py @@ -21,9 +21,10 @@ import time from scipy.signal import convolve MAX_DISTANCE = 200.0 -FRICTION_COEFFICENT = 0.9 +FRICTION_COEFFICENT = 1 EPSILON = 0.00001 + def value_to_color(v, min_value, max_value): r = g = b = 0.0 scope = max_value - min_value @@ -32,12 +33,13 @@ def value_to_color(v, min_value, max_value): b = max(0, 1.0 - abs(2.0 * normalized)) g = max(0, 1.0 - abs(2.0 * normalized - 1.0)) r = max(0, 1.0 - abs(2.0 * normalized - 2.0)) - l = np.sqrt((r*r + b*b + g*g)) + l = np.sqrt((r * r + b * b + g * g)) r /= l g /= l b /= l return r, g, b + class Client: def __init__(self, test=False, pos=[0, 0, 0]): with open('./config.json', 'r') as f: @@ -90,7 +92,8 @@ class Client: glAttachShader(self.normal_program[key], self.fragment_shader_id) 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) for x_pos in range(0, 100): @@ -113,11 +116,44 @@ class Client: 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.v_map_y = np.zeros((100, 100, 1)) - self.v_map_z = np.zeros((100, 100, 1)) + self.relaxation_time = 0.55 # 0.55 + self.w_a = [ + 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: glutMainLoop() @@ -125,8 +161,6 @@ class Client: self.display() self.resize(100, 100) - - def display(self): glClearColor(0, 0, 0, 0) glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT) @@ -172,171 +206,143 @@ class Client: glUniform1f(widthid, self.width) 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() glutSwapBuffers() - max_value = np.max(self.heat_map) - # min_value = np.min(self.heat_map) 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)) - max_value = np.max(vel) - min_value = np.min(vel) + vel = np.sqrt(np.sum(np.square(self.u), axis=3)) + max_value_vel = np.max(vel) + # max_value_vel = np.sqrt(3) + + print('round') + print(max_value_n) for x_pos in range(0, 100): for y_pos in range(0, 100): 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(vel[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_vel) self.world_provider.world.set_color(x_pos, y_pos, z_pos, r, g, b) - # friction - # self.heat_map += np.sqrt(np.square(self.v_map_x * (1.0 - FRICTION_COEFFICENT)) + - # np.square(self.v_map_y * (1.0 - FRICTION_COEFFICENT)) + - # 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])))) + self.world_provider.world.set_color(int(round(self.test_pixel[0])), + int(round(self.test_pixel[1])), + int(round(self.test_pixel[2])), 1.0, 1.0, 1.0) - friction_dispersion = (1.0 -FRICTION_COEFFICENT) / 4 - # 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 + # self.u *= 0.95 - if target_z + 1 < self.v_map_x.shape[2] - 1: - new_v_map_z[target_x, target_y, target_z + 1] += 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 > 0: - new_v_map_z[target_x, target_y, target_z - 1] -= 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 - # y - new_v_map_y[target_x, target_y, target_z] += self.v_map_y[x_pos, y_pos, z_pos] * FRICTION_COEFFICENT - 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_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 + old_n_sum = np.sum(self.n) + for a in range(len(self.e_a)): + e_au = np.sum(self.e_a[a] * self.u, axis=3) + uu = np.sum(self.u * self.u, axis=3) + 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) + print(np.max(self.n_a_eq[0])) + if not self.compressible: + excess = (self.n_a_eq[0] > self.max_n) * (self.n_a_eq[0] - self.max_n) + dir_sum = np.sum(self.n_a_eq[1:], axis=0) + self.n_a_eq[1:] += excess * 1/8 #(self.n_a_eq[1:] / dir_sum) + self.n_a_eq[0] -= excess - if target_z + 1 < self.v_map_x.shape[2] - 1: - new_v_map_z[target_x, target_y, target_z + 1] += 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 > 0: - new_v_map_z[target_x, target_y, target_z - 1] -= 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 - # z - new_v_map_z[target_x, target_y, target_z] += self.v_map_z[x_pos, y_pos, z_pos] * FRICTION_COEFFICENT - 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 - 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_x - 1 > 0: - new_v_map_x[target_x - 1, target_y, target_z] -= self.v_map_z[x_pos, y_pos, z_pos] * friction_dispersion - else: - new_v_map_z[target_x, target_y, target_z] += self.v_map_z[x_pos, y_pos, z_pos] * friction_dispersion + n_a_t_1 = np.zeros((len(self.e_a),) + self.field) + 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] - 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_z[x_pos, y_pos, z_pos] * friction_dispersion - 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 > 0: - new_v_map_y[target_x, target_y - 1, target_z] -= self.v_map_z[x_pos, y_pos, z_pos] * friction_dispersion - else: - 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 + clipped_dir = np.zeros(3) + clipped_dir[index] = self.e_a[a][index] - filter_mat = np.array([[-1.0], [0], [1.0]]) / 2.0 - new_v_map_x[:, 0, :] += convolve(new_v_map_y[:, 0, :], filter_mat, 'same') - 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 + e_a_index = np.where(np.all(self.e_a == e_a_clipped, axis=1))[0][0] + # e_a_index = 0 - # corners - new_v_map_x[0, 0, 0] = new_v_map_y[0, 0, 0] = new_v_map_z[0, 0, 0] = 0 - new_v_map_x[-1, 0, 0] = new_v_map_y[-1, 0, 0] = new_v_map_z[-1, 0, 0] = 0 - new_v_map_x[-1, -1, 0] = new_v_map_y[-1, -1, 0] = new_v_map_z[-1, -1, 0] = 0 - new_v_map_x[-1, -1, -1] = new_v_map_y[-1, -1, -1] = new_v_map_z[-1, -1, -1] = 0 - new_v_map_x[0, -1, -1] = new_v_map_y[0, -1, -1] = new_v_map_z[0, -1, -1] = 0 - new_v_map_x[0, -1, 0] = new_v_map_y[0, -1, 0] = new_v_map_z[0, -1, 0] = 0 - new_v_map_x[-1, -1, 0] = new_v_map_y[-1, -1, 0] = new_v_map_z[-1, -1, 0] = 0 - new_v_map_x[-1, 0, -1] = new_v_map_y[-1, 0, -1] = new_v_map_z[-1, 0, -1] = 0 + if index == 0: + if self.e_a[a][index] > 0: + slice_index = -1 + else: + slice_index = 0 + n_a_t_1[e_a_index, slice_index, :, :] += temp[slice_index, :, :] + temp[slice_index, :, :] *= 0 - self.v_map_x = new_v_map_x - self.v_map_y = new_v_map_y - self.v_map_z = new_v_map_z + if index == 1: + if self.e_a[a][index] > 0: + slice_index = -1 + else: + slice_index = 0 + n_a_t_1[e_a_index, :, slice_index, :] += temp[:, slice_index, :] + temp[:, slice_index, :] *= 0 - filter_mat = (np.zeros((3, 3, 1)) + 1.0) / 9.0 + if index == 2: + if self.e_a[a][index] > 0: + slice_index = -1 + else: + slice_index = 0 + n_a_t_1[e_a_index, :, :, slice_index] += temp[:, :, slice_index] + temp[:, :, slice_index] *= 0 - v_map = np.pad(self.v_map_x, 1, 'edge') - self.v_map_x = convolve(v_map, filter_mat, mode='same')[1:-1, 1:-1, 1:2] - v_map = np.pad(self.v_map_y, 1, 'edge') - self.v_map_y = convolve(v_map, filter_mat, mode='same')[1:-1, 1:-1, 1:2] - v_map = np.pad(self.v_map_z, 1, 'edge') - 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.n_a = n_a_t_1 + if np.min(self.n_a) < 0: + test = 1 + self.n_a = np.maximum(0, self.n_a) + self.n_a = old_n_sum * self.n_a / np.sum(self.n_a) + self.n = np.sum(self.n_a, axis=0, keepdims=False) + # self.n = np.sum(np.abs(n_a_t_1), axis=0, keepdims=False) - self.heat_map = heat_map - # dispersing heat - heat_map = np.pad(self.heat_map, 1, 'edge') - self.heat_map = convolve(heat_map, filter_mat, mode='same')[1:-1, 1:-1, 1:2] - # 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) + # stabilise the number because of rounding errors + self.n = old_n_sum * self.n / np.sum(self.n) - print(1.0 / (time.time() - self.time)) + self.u *= 0 + for a in range(len(self.e_a)): + self.u[:, :, :, 0] += self.n_a[a] * self.e_a[a][0] + self.u[:, :, :, 1] += self.n_a[a] * self.e_a[a][1] + self.u[:, :, :, 2] += self.n_a[a] * self.e_a[a][2] + self.u[:, :, :, 0] /= self.n + self.u[:, :, :, 1] /= self.n + self.u[:, :, :, 2] /= self.n + + self.u[self.n == 0] = 0 + + length = np.sqrt(np.sum(np.square(self.u), axis=3, keepdims=True)) + + # gravity + gravity_applies = self.n < self.w_a[0] + gravity_applies[:, :-1, :] = gravity_applies[:, 1:, :] + gravity_applies[:, -1, :] = False + self.u[gravity_applies, 1] -= 0.01 + + new_lengths = np.sqrt(np.sum(np.square(self.u), axis=3, keepdims=True)) + self.u = self.u / new_lengths * length + zero_length = (new_lengths == 0) + self.u[zero_length[:, :, :, 0], :] = 0 + + 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() glutPostRedisplay() @@ -382,5 +388,6 @@ class Client: glutFullScreenToggle() # print(key) + if __name__ == '__main__': client = Client(pos=[-50, -50, -200])