Skip to content
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
48 changes: 48 additions & 0 deletions examples/RigidBallSoftTrampoline/collision_force.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,48 @@
import numpy as np
import scipy.sparse as spp


def collision_force(k, phi_1, phi_2, blend_w, grid_size_z, grid_size_r, eps, dx):
"""
computes the collision force between 2 bodies
using their level sets
"""

gradient = (
spp.diags(
[-0.5, 0.5, -0.5, 0.5],
[-1, 1, grid_size_r - 1, -grid_size_r + 1],
shape=(grid_size_r, grid_size_r),
format="csr",
)
/ dx
)

gradient_T = (
spp.diags(
[-0.5, 0.5, -0.5, 0.5],
[-1, 1, grid_size_z - 1, -grid_size_z + 1],
shape=(grid_size_z, grid_size_z),
format="csr",
)
/ dx
)

phi_diff = 0.5 * (phi_2 - phi_1)
delta_phi = phi_diff * 0
delta_phi += (
(np.fabs(phi_diff) < blend_w)
* 0.5
* (1 + np.cos(np.pi * phi_diff / blend_w))
/ blend_w
)
phi_diff_x = np.mat(phi_diff) * gradient_T
phi_diff_y = gradient * np.mat(phi_diff)
s = np.sqrt(np.square(phi_diff_x) + np.square(phi_diff_y))
phi_diff_x /= s + eps
phi_diff_y /= s + eps
S = phi_diff / np.sqrt(phi_diff**2 + eps**2)
body_mask = np.logical_or((phi_1 > 0), (phi_2 > 0))
force_x = body_mask * S * phi_diff_x * k * delta_phi
force_y = body_mask * S * phi_diff_y * k * delta_phi
return force_x, force_y
Comment on lines +5 to +48
Copy link
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

please write the vectorized versions (slice-based, similar to functions in kernels), Let's skip using the sparse matrices if possible. Also, collision forces should be a part of the main package.

31 changes: 31 additions & 0 deletions examples/RigidBallSoftTrampoline/curl.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,31 @@
import numpy as np
import scipy.sparse as spp


def curl(ax, ay, grid_size_z, grid_size_r, dx):
"""
compute curl of the vector using CD
"""

gradient = (
spp.diags(
[-0.5, 0.5, -0.5, 0.5],
[-1, 1, grid_size_r - 1, -grid_size_r + 1],
shape=(grid_size_r, grid_size_r),
format="csr",
)
/ dx
)

gradient_T = (
spp.diags(
[-0.5, 0.5, -0.5, 0.5],
[-1, 1, grid_size_z - 1, -grid_size_z + 1],
shape=(grid_size_z, grid_size_z),
format="csr",
)
/ dx
)

vort = np.mat(ay) * gradient_T - gradient * np.mat(ax)
Comment on lines +5 to +30
Copy link
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

can't we use compute_vorticity_from_velocity for this? maybe a good idea is to rename that function as compute_vorticity_from_velocity_forcing?

return vort
Loading