From 9cda7396983328df0b9ffc5e24cfd155361653a1 Mon Sep 17 00:00:00 2001 From: zomseffen Date: Sat, 13 Feb 2021 21:31:39 +0100 Subject: [PATCH] fluid sim --- Client/Client.py | 184 +++++++++++++++++++++++++++++++++++++++++++++-- 1 file changed, 178 insertions(+), 6 deletions(-) diff --git a/Client/Client.py b/Client/Client.py index 9baea65..0aa5dd4 100644 --- a/Client/Client.py +++ b/Client/Client.py @@ -119,7 +119,7 @@ class Client: self.field = (100, 100, 1) self.e_a = np.array([ - [0, 0, 0], + # [0, 0, 0], [1, 0, 0], [1, 1, 0], [0, 1, 0], @@ -145,12 +145,20 @@ class Client: 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.zeros(self.field, np.int) + self.n[:, 50:, :] += 1 + self.viscosity = np.reshape(self.n * 0.01, self.field + (1,)) self.gravity_applies = np.zeros(self.field) # 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.f = np.zeros(self.field + (self.e_a.shape[1],)) + self.true_pos = np.zeros(self.field + (self.e_a.shape[1],)) + self.pos_indices = np.zeros(self.field + (self.e_a.shape[1],), dtype=np.int) + for x in range(self.field[0]): + for y in range(self.field[1]): + for z in range(self.field[2]): + self.pos_indices[x, y, z, :] = [x, y, z] self.compressible = True self.max_n = self.w_a[0] @@ -218,7 +226,7 @@ class Client: max_value_n = np.max(self.n) # max_value_n = 1.0 - vel = np.sqrt(np.sum(np.square(self.u), axis=3)) *self.n + vel = np.sqrt(np.sum(np.square(self.u), axis=3)) * self.n max_value_vel = np.max(vel) # max_value_vel = np.sqrt(3) @@ -245,7 +253,171 @@ class Client: int(round(self.test_pixel[1])), int(round(self.test_pixel[2])), 1.0, 1.0, 1.0) - # print(1.0 / (time.time() - self.time)) + new_f = np.zeros(self.f.shape) + # viscosity + neighbours_filter = np.zeros((3, 3, 3), np.int) + 1 + neighbours = convolve(np.pad(self.n, 1, constant_values=0), neighbours_filter, mode='valid') * self.n + neighbours = np.reshape(neighbours, self.field + (1,)) + forces_to_share_per_neighbor = self.f * self.viscosity + length = np.sqrt(np.sum(np.square(forces_to_share_per_neighbor), axis=3, keepdims=True)) + direction = forces_to_share_per_neighbor / length + direction[(length == 0)[:, :, :, 0]] = 0 + ############################## experimental + for a in range(len(self.e_a)): + unit = self.e_a[a] / np.sqrt(np.sum(np.square(self.e_a[a]))) + scalar = np.sum(direction * unit, axis=3, keepdims=True) + altered_direction = direction + scalar * unit + altered_length = np.sqrt(np.sum(np.square(altered_direction), axis=3, keepdims=True)) + altered_direction = altered_direction / altered_length + altered_direction[(altered_length == 0)[:, :, :, 0]] = 0 + + f_2_add = length * altered_direction + + new_f[max(0, self.e_a[a][0]):min(new_f.shape[0], new_f.shape[0] + self.e_a[a][0]), + max(0, self.e_a[a][1]):min(new_f.shape[1], new_f.shape[1] + self.e_a[a][1]), + max(0, self.e_a[a][2]):min(new_f.shape[2], new_f.shape[2] + self.e_a[a][2])] += f_2_add[ + max(0, -self.e_a[a][0]):min( + new_f.shape[0], + new_f.shape[0] - + self.e_a[a][0]), + max(0, -self.e_a[a][1]):min( + new_f.shape[1], + new_f.shape[1] - + self.e_a[a][1]), + max(0, -self.e_a[a][2]):min( + new_f.shape[2], + new_f.shape[2] - + self.e_a[a][2])] + + new_f += self.f - forces_to_share_per_neighbor * neighbours + ######################################### + new_f += self.f + # TODO movement generating things + # gravity + new_f[:, :, :, 1] -= 1.0 * self.n + # clean up + new_f[self.n == 0] = 0 + + self.f = new_f + # collision and movement + collision = True + iterations = 0 + while collision: + f_length = np.sqrt(np.sum(np.square(self.f), axis=3, keepdims=True)) + f_direction = self.f / f_length + f_direction[f_length[:, :, :, 0] == 0] = 0 + velocity = f_direction * np.sqrt(f_length) / np.reshape(self.n, self.field + (1,)) # TODO replace self.n by mass + velocity[self.n == 0] = 0 + timestep = min(1, 0.5 / np.max(np.sqrt(np.sum(np.square(velocity), axis=3)))) + if iterations > 20: + print('Takes too long!') + timestep /= 10 + new_pos = self.true_pos + velocity * timestep + new_pos_round = np.round(new_pos) + moved = np.logical_or.reduce(new_pos_round != [0, 0, 0], axis=3) + pos_change_targets = new_pos_round[moved] + self.pos_indices[moved] + + # handle bordercases + bordercase = np.array(moved) + bordercase[moved] = np.logical_or( + np.logical_or(np.logical_or(0 > pos_change_targets[:, 0], pos_change_targets[:, 0] >= self.field[0]), + np.logical_or(0 > pos_change_targets[:, 1], pos_change_targets[:, 1] >= self.field[1])), + np.logical_or(0 > pos_change_targets[:, 2], pos_change_targets[:, 2] >= self.field[2])) + self.f[bordercase] *= -1 + velocity[bordercase] *= -1 + + # recalculate targets + new_pos = self.true_pos + velocity * timestep + new_pos_round = np.round(new_pos) + new_pos_target = new_pos_round + self.pos_indices + + contenders = np.zeros(self.field) + starts = self.pos_indices[self.n != 0] + targets = new_pos_target[self.n != 0].astype(np.int) + speeds = velocity[self.n != 0] + forces = new_f[self.n != 0] + max_speeds = np.zeros(self.field) + fast_pos = {} + is_stayer = [] + for index in range(len(targets)): + target = targets[index] + speed = np.sqrt(np.sum(np.square(speeds[index]))) + contenders[target[0], target[1], target[2]] += 1 + + # new max speed and we do not update stayers + if speed > max_speeds[target[0], target[1], target[2]] and tuple(target) not in is_stayer: + max_speeds[target[0], target[1], target[2]] = speed + fast_pos[tuple(target)] = index + + # atoms that are already there are there (and stay there) are the fastest + start = starts[index] + if np.all(start == target): + fast_pos[tuple(target)] = index + is_stayer.append(tuple(target)) + + collision = np.max(contenders) > 1 + + # we only need collision if there is one + if collision: + # go through the movers again + for index in range(len(targets)): + target = targets[index] + # collision here? + if contenders[target[0], target[1], target[2]] > 1: + # the fastest one does not need to do anything + if index != fast_pos[tuple(target)]: + force = forces[index] + fastest = fast_pos[tuple(target)] + # TODO use relative weight? + forces[fastest] += 0.5*force + forces[index] *= 0.5 + new_f[self.n != 0] = forces + + + self.f = new_f + + iterations += 1 + + # final calculation + f_length = np.sqrt(np.sum(np.square(self.f), axis=3, keepdims=True)) + f_direction = self.f / f_length + f_direction[f_length[:, :, :, 0] == 0] = 0 + velocity = f_direction * np.sqrt(f_length) / np.reshape(self.n, self.field + (1,)) # TODO replace self.n by mass + velocity[self.n == 0] = 0 + #timestep = min(1, 0.5 / np.max(np.sqrt(np.sum(np.square(velocity), axis=3)))) + new_pos = self.true_pos + velocity * timestep + new_pos_round = np.round(new_pos) + moved = np.logical_or.reduce(new_pos_round[:, :, :] != [0, 0, 0], axis=3) + not_moved = np.logical_and.reduce(new_pos_round[:, :, :] == [0, 0, 0], axis=3) + pos_change_targets = (new_pos_round[moved] + self.pos_indices[moved]).astype(np.int) + movers = self.pos_indices[moved] + + print('timestep: %f' % timestep) + + update_n = np.zeros(self.n.shape, np.int) + update_f = np.zeros(self.f.shape) + update_u = np.zeros(self.u.shape) + update_true_pos = np.zeros(self.true_pos.shape) + + update_n[not_moved] = self.n[not_moved] + update_f[not_moved, :] = self.f[not_moved, :] + update_u[not_moved, :] = velocity[not_moved, :] + update_true_pos[not_moved, :] = new_pos[not_moved, :] + + for indice in range(len(movers)): + mover = movers[indice] + target = pos_change_targets[indice] + update_n[target[0], target[1], target[2]] = self.n[mover[0], mover[1], mover[2]] + update_f[target[0], target[1], target[2], :] = self.f[mover[0], mover[1], mover[2], :] + update_u[target[0], target[1], target[2], :] = velocity[mover[0], mover[1], mover[2], :] + update_true_pos[target[0], target[1], target[2], :] = new_pos[mover[0], mover[1], mover[2], :] - new_pos_round[mover[0], mover[1], mover[2], :] + + self.n = update_n + self.f = update_f + self.u = update_u + self.true_pos = update_true_pos + + print(1.0 / (time.time() - self.time)) self.time = time.time() glutPostRedisplay() @@ -282,7 +454,7 @@ class Client: self.opening += 0.25 if key == b'+': - self.state = (self.state +1) % 3 + self.state = (self.state + 1) % 3 if key == b'r': print(self.cx, self.cy, self.opening)