Skip to content

Commit cba2f87

Browse files
liminchenliminchen
committed
moving DBC
1 parent ae2d997 commit cba2f87

3 files changed

Lines changed: 46 additions & 29 deletions

File tree

‎mov_dirichlet/SpringEnergy.py‎

Lines changed: 3 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -1,21 +1,19 @@
11
import numpy as np
22

3-
k = 10
4-
5-
def val(x, m, DBC, DBC_target):
3+
def val(x, m, DBC, DBC_target, k):
64
sum = 0.0
75
for i in range(0, len(DBC)):
86
diff = x[DBC[i]] - DBC_target[i]
97
sum += 0.5 * k * m[DBC[i]] * diff.dot(diff)
108
return sum
119

12-
def grad(x, m, DBC, DBC_target):
10+
def grad(x, m, DBC, DBC_target, k):
1311
g = np.array([[0.0, 0.0]] * len(x))
1412
for i in range(0, len(DBC)):
1513
g[DBC[i]] = k * m[DBC[i]] * (x[DBC[i]] - DBC_target[i])
1614
return g
1715

18-
def hess(x, m, DBC, DBC_target):
16+
def hess(x, m, DBC, DBC_target, k):
1917
IJV = [[0] * 0, [0] * 0, np.array([0.0] * 0)]
2018
for i in range(0, len(DBC)):
2119
for d in range(0, 2):

‎mov_dirichlet/simulator.py‎

Lines changed: 6 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -13,10 +13,12 @@
1313
k = 2e4 # spring stiffness
1414
n_seg = 4 # num of segments per side of the square
1515
h = 0.01 # time step size in s
16-
DBC = [] # no nodes need to be fixed
17-
ground_n = np.array([0.0, 1.0]) # normal of the slope
16+
DBC = [(n_seg + 1) * (n_seg + 1)] # dirichlet node index
17+
DBC_v = [np.array([0.0, -0.5])] # dirichlet node velocity
18+
DBC_limit = [np.array([0.0, -0.6])] # dirichlet node limit position
19+
ground_n = np.array([0.0, 1.0]) # normal of the slope
1820
ground_n /= np.linalg.norm(ground_n) # normalize ground normal vector just in case
19-
ground_o = np.array([0.0, -1.0]) # a point on the slope
21+
ground_o = np.array([0.0, -1.0]) # a point on the slope
2022
mu = 0.11 # friction coefficient of the slope
2123

2224
# initialize simulation
@@ -69,7 +71,7 @@ def screen_projection(x):
6971
pygame.display.flip() # flip the display
7072

7173
# step forward simulation and wait for screen refresh
72-
[x, v] = time_integrator.step_forward(x, e, v, m, l2, k, ground_n, ground_o, contact_area, mu, is_DBC, h, 1e-2)
74+
[x, v] = time_integrator.step_forward(x, e, v, m, l2, k, ground_n, ground_o, contact_area, mu, is_DBC, DBC, DBC_v, DBC_limit, h, 1e-2)
7375
time_step += 1
7476
pygame.time.wait(int(h * 1000))
7577

‎mov_dirichlet/time_integrator.py‎

Lines changed: 37 additions & 20 deletions
Original file line numberDiff line numberDiff line change
@@ -13,55 +13,67 @@
1313
import FrictionEnergy
1414
import SpringEnergy
1515

16-
def step_forward(x, e, v, m, l2, k, n, o, contact_area, mu, is_DBC, h, tol):
16+
def step_forward(x, e, v, m, l2, k, n, o, contact_area, mu, is_DBC, DBC, DBC_v, DBC_limit, h, tol):
1717
x_tilde = x + v * h # implicit Euler predictive position
1818
x_n = copy.deepcopy(x)
1919
mu_lambda = BarrierEnergy.compute_mu_lambda(x, n, o, contact_area, mu) # compute mu * lambda for each node using x^n
20+
DBC_target = [] # target position of each DBC in the current time step
21+
for i in range(0, len(DBC)):
22+
if (DBC_limit[i] - x_n[DBC[i]]).dot(DBC_v[i]) > 0:
23+
DBC_target.append(x_n[DBC[i]] + h * DBC_v[i])
24+
else:
25+
DBC_target.append(x_n[DBC[i]])
26+
DBC_stiff = 10 # initialize stiffness for DBC springs
2027

2128
# Newton loop
2229
iter = 0
23-
E_last = IP_val(x, e, x_tilde, m, l2, k, n, o, contact_area, (x - x_n) / h, mu_lambda, h)
24-
p = search_dir(x, e, x_tilde, m, l2, k, n, o, contact_area, (x - x_n) / h, mu_lambda, is_DBC, h)
25-
while LA.norm(p, inf) / h > tol:
30+
E_last = IP_val(x, e, x_tilde, m, l2, k, n, o, contact_area, (x - x_n) / h, mu_lambda, DBC, DBC_target, DBC_stiff, h)
31+
[p, DBC_satisfied] = search_dir(x, e, x_tilde, m, l2, k, n, o, contact_area, (x - x_n) / h, mu_lambda, is_DBC, DBC, DBC_target, DBC_stiff, tol, h)
32+
while (LA.norm(p, inf) / h > tol) | (sum(DBC_satisfied) != len(DBC)): # also check whether all DBCs are satisfied
2633
print('Iteration', iter, ':')
2734
print('residual =', LA.norm(p, inf) / h)
2835

36+
if (LA.norm(p, inf) / h <= tol) & (sum(DBC_satisfied) != len(DBC)):
37+
# increase DBC stiffness and recompute energy value record
38+
DBC_stiff *= 2
39+
E_last = IP_val(x, e, x_tilde, m, l2, k, n, o, contact_area, (x - x_n) / h, mu_lambda, DBC, DBC_target, DBC_stiff, h)
40+
2941
# filter line search
3042
alpha = BarrierEnergy.init_step_size(x, n, o, p) # avoid interpenetration and tunneling
31-
while IP_val(x + alpha * p, e, x_tilde, m, l2, k, n, o, contact_area, (x - x_n) / h, mu_lambda, h) > E_last:
43+
while IP_val(x + alpha * p, e, x_tilde, m, l2, k, n, o, contact_area, (x - x_n) / h, mu_lambda, DBC, DBC_target, DBC_stiff, h) > E_last:
3244
alpha /= 2
3345
print('step size =', alpha)
3446

3547
x += alpha * p
36-
E_last = IP_val(x, e, x_tilde, m, l2, k, n, o, contact_area, (x - x_n) / h, mu_lambda, h)
37-
p = search_dir(x, e, x_tilde, m, l2, k, n, o, contact_area, (x - x_n) / h, mu_lambda, is_DBC, h)
48+
E_last = IP_val(x, e, x_tilde, m, l2, k, n, o, contact_area, (x - x_n) / h, mu_lambda, DBC, DBC_target, DBC_stiff, h)
49+
[p, DBC_satisfied] = search_dir(x, e, x_tilde, m, l2, k, n, o, contact_area, (x - x_n) / h, mu_lambda, is_DBC, DBC, DBC_target, DBC_stiff, tol, h)
3850
iter += 1
3951

4052
v = (x - x_n) / h # implicit Euler velocity update
4153
return [x, v]
4254

43-
def IP_val(x, e, x_tilde, m, l2, k, n, o, contact_area, v, mu_lambda, h):
55+
def IP_val(x, e, x_tilde, m, l2, k, n, o, contact_area, v, mu_lambda, DBC, DBC_target, DBC_stiff, h):
4456
return InertiaEnergy.val(x, x_tilde, m) + h * h * ( # implicit Euler
4557
MassSpringEnergy.val(x, e, l2, k) +
4658
GravityEnergy.val(x, m) +
4759
BarrierEnergy.val(x, n, o, contact_area) +
4860
FrictionEnergy.val(v, mu_lambda, h, n)
49-
) + SpringEnergy.val(x, m, [len(x) - 1], [np.array([0.0, 0.6])])
61+
) + SpringEnergy.val(x, m, DBC, DBC_target, DBC_stiff)
5062

51-
def IP_grad(x, e, x_tilde, m, l2, k, n, o, contact_area, v, mu_lambda, h):
63+
def IP_grad(x, e, x_tilde, m, l2, k, n, o, contact_area, v, mu_lambda, DBC, DBC_target, DBC_stiff, h):
5264
return InertiaEnergy.grad(x, x_tilde, m) + h * h * ( # implicit Euler
5365
MassSpringEnergy.grad(x, e, l2, k) +
5466
GravityEnergy.grad(x, m) +
5567
BarrierEnergy.grad(x, n, o, contact_area) +
5668
FrictionEnergy.grad(v, mu_lambda, h, n)
57-
) + SpringEnergy.grad(x, m, [len(x) - 1], [np.array([0.0, 0.6])])
69+
) + SpringEnergy.grad(x, m, DBC, DBC_target, DBC_stiff)
5870

59-
def IP_hess(x, e, x_tilde, m, l2, k, n, o, contact_area, v, mu_lambda, h):
71+
def IP_hess(x, e, x_tilde, m, l2, k, n, o, contact_area, v, mu_lambda, DBC, DBC_target, DBC_stiff, h):
6072
IJV_In = InertiaEnergy.hess(x, x_tilde, m)
6173
IJV_MS = MassSpringEnergy.hess(x, e, l2, k)
6274
IJV_B = BarrierEnergy.hess(x, n, o, contact_area)
6375
IJV_F = FrictionEnergy.hess(v, mu_lambda, h, n)
64-
IJV_S = SpringEnergy.hess(x, m, [len(x) - 1], [np.array([0.0, 0.6])])
76+
IJV_S = SpringEnergy.hess(x, m, DBC, DBC_target, DBC_stiff)
6577
IJV_MS[2] *= h * h # implicit Euler
6678
IJV_B[2] *= h * h # implicit Euler
6779
IJV_F[2] *= h * h # implicit Euler
@@ -72,14 +84,19 @@ def IP_hess(x, e, x_tilde, m, l2, k, n, o, contact_area, v, mu_lambda, h):
7284
H = sparse.coo_matrix((IJV[2], (IJV[0], IJV[1])), shape=(len(x) * 2, len(x) * 2)).tocsr()
7385
return H
7486

75-
def search_dir(x, e, x_tilde, m, l2, k, n, o, contact_area, v, mu_lambda, is_DBC, h):
76-
projected_hess = IP_hess(x, e, x_tilde, m, l2, k, n, o, contact_area, v, mu_lambda, h)
77-
reshaped_grad = IP_grad(x, e, x_tilde, m, l2, k, n, o, contact_area, v, mu_lambda, h).reshape(len(x) * 2, 1)
78-
# eliminate DOF by modifying gradient and Hessian for DBC:
87+
def search_dir(x, e, x_tilde, m, l2, k, n, o, contact_area, v, mu_lambda, is_DBC, DBC, DBC_target, DBC_stiff, tol, h):
88+
projected_hess = IP_hess(x, e, x_tilde, m, l2, k, n, o, contact_area, v, mu_lambda, DBC, DBC_target, DBC_stiff, h)
89+
reshaped_grad = IP_grad(x, e, x_tilde, m, l2, k, n, o, contact_area, v, mu_lambda, DBC, DBC_target, DBC_stiff, h).reshape(len(x) * 2, 1)
90+
# check whether each DBC is satisfied
91+
DBC_satisfied = [False] * len(x)
92+
for i in range(0, len(DBC)):
93+
if LA.norm(x[DBC[i]] - DBC_target[i]) / h < tol:
94+
DBC_satisfied[DBC[i]] = True
95+
# eliminate DOF if it's a satisfied DBC by modifying gradient and Hessian for DBC:
7996
for i, j in zip(*projected_hess.nonzero()):
80-
if is_DBC[int(i / 2)] | is_DBC[int(j / 2)]:
97+
if (is_DBC[int(i / 2)] & DBC_satisfied[int(i / 2)]) | (is_DBC[int(j / 2)] & DBC_satisfied[int(i / 2)]):
8198
projected_hess[i, j] = (i == j)
8299
for i in range(0, len(x)):
83-
if is_DBC[i]:
100+
if is_DBC[i] & DBC_satisfied[i]:
84101
reshaped_grad[i * 2] = reshaped_grad[i * 2 + 1] = 0.0
85-
return spsolve(projected_hess, -reshaped_grad).reshape(len(x), 2)
102+
return [spsolve(projected_hess, -reshaped_grad).reshape(len(x), 2), DBC_satisfied]

0 commit comments

Comments
 (0)