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

fluid sim

This commit is contained in:
zomseffen 2021-02-13 21:31:39 +01:00
parent afc0fff4fa
commit 9cda739698
1 changed files with 178 additions and 6 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

184
Client/Client.py
View File

@ -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)