Skip to content

Commit 0b92dff

Browse files
committed
added anchors
1 parent 2a93cfe commit 0b92dff

2 files changed

Lines changed: 9 additions & 1 deletion

File tree

2_dirichlet/simulator.py

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -13,7 +13,9 @@
1313
k = 1e3 # spring stiffness
1414
n_seg = 4 # num of segments per side of the square
1515
h = 0.02 # time step size in s
16+
# ANCHOR: DBC_def
1617
DBC = [n_seg, (n_seg + 1) * (n_seg + 1) - 1] # fix the left and right top nodes
18+
# ANCHOR_END: DBC_def
1719

1820
# initialize simulation
1921
[x, e] = square_mesh.generate(side_len, n_seg) # node positions and edge node indices
@@ -25,10 +27,12 @@
2527
diff = x[e[i][0]] - x[e[i][1]]
2628
l2.append(diff.dot(diff))
2729
k = [k] * len(e) # spring stiffness
30+
# ANCHOR: DBC_mask
2831
# identify whether a node is Dirichlet
2932
is_DBC = [False] * len(x)
3033
for i in DBC:
3134
is_DBC[i] = True
35+
# ANCHOR_END: DBC_mask
3236
# simulation with visualization
3337
resolution = np.array([900, 900])
3438
offset = resolution / 2

2_dirichlet/time_integrator.py

Lines changed: 5 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -36,11 +36,13 @@ def step_forward(x, e, v, m, l2, k, is_DBC, h, tol):
3636
v = (x - x_n) / h # implicit Euler velocity update
3737
return [x, v]
3838

39+
# ANCHOR: ADDING_GRAVITY
3940
def IP_val(x, e, x_tilde, m, l2, k, h):
4041
return InertiaEnergy.val(x, x_tilde, m) + h * h * (MassSpringEnergy.val(x, e, l2, k) + GravityEnergy.val(x, m)) # implicit Euler
4142

4243
def IP_grad(x, e, x_tilde, m, l2, k, h):
4344
return InertiaEnergy.grad(x, x_tilde, m) + h * h * (MassSpringEnergy.grad(x, e, l2, k) + GravityEnergy.grad(x, m)) # implicit Euler
45+
# ANCHOR_END: ADDING_GRAVITY
4446

4547
def IP_hess(x, e, x_tilde, m, l2, k, h):
4648
IJV_In = InertiaEnergy.hess(x, x_tilde, m)
@@ -50,6 +52,7 @@ def IP_hess(x, e, x_tilde, m, l2, k, h):
5052
H = sparse.coo_matrix((IJV[2], (IJV[0], IJV[1])), shape=(len(x) * 2, len(x) * 2)).tocsr()
5153
return H
5254

55+
# ANCHOR: search_dir
5356
def search_dir(x, e, x_tilde, m, l2, k, is_DBC, h):
5457
projected_hess = IP_hess(x, e, x_tilde, m, l2, k, h)
5558
reshaped_grad = IP_grad(x, e, x_tilde, m, l2, k, h).reshape(len(x) * 2, 1)
@@ -60,4 +63,5 @@ def search_dir(x, e, x_tilde, m, l2, k, is_DBC, h):
6063
for i in range(0, len(x)):
6164
if is_DBC[i]:
6265
reshaped_grad[i * 2] = reshaped_grad[i * 2 + 1] = 0.0
63-
return spsolve(projected_hess, -reshaped_grad).reshape(len(x), 2)
66+
return spsolve(projected_hess, -reshaped_grad).reshape(len(x), 2)
67+
#ANCHOR_END: search_dir

0 commit comments

Comments
 (0)