You signed in with another tab or window. Reload to refresh your session.You signed out in another tab or window. Reload to refresh your session.You switched accounts on another tab or window. Reload to refresh your session.Dismiss alert
Copy file name to clipboardExpand all lines: 6_inv_free/time_integrator.py
+8-10Lines changed: 8 additions & 10 deletions
Original file line number
Diff line number
Diff line change
@@ -16,7 +16,7 @@
16
16
17
17
18
18
19
-
defstep_forward(x, e, v, m, vol, IB, mu_lame, lam, n, o, contact_area, mu, is_DBC, DBC, DBC_v, DBC_limit, DBC_stiff_,h, tol):
19
+
defstep_forward(x, e, v, m, vol, IB, mu_lame, lam, n, o, contact_area, mu, is_DBC, DBC, DBC_v, DBC_limit, DBC_stiff, h, tol):
20
20
x_tilde=x+v*h# implicit Euler predictive position
21
21
x_n=copy.deepcopy(x)
22
22
mu_lambda=BarrierEnergy.compute_mu_lambda(x, n, o, contact_area, mu) # compute mu * lambda for each node using x^n
@@ -26,33 +26,31 @@ def step_forward(x, e, v, m, vol, IB, mu_lame, lam, n, o, contact_area, mu, is_D
26
26
DBC_target.append(x_n[DBC[i]] +h*DBC_v[i])
27
27
else:
28
28
DBC_target.append(x_n[DBC[i]])
29
-
DBC_stiff=DBC_stiff_[0]
30
29
31
30
# Newton loop
32
31
iter=0
33
-
E_last=IP_val(x, e, x_tilde, m, vol, IB, mu_lame, lam, n, o, contact_area, (x-x_n) /h, mu_lambda, DBC, DBC_target, DBC_stiff, h)
34
-
[p, DBC_satisfied] =search_dir(x, e, x_tilde, m, vol, IB, mu_lame, lam, n, o, contact_area, (x-x_n) /h, mu_lambda, is_DBC, DBC, DBC_target, DBC_stiff, tol, h)
32
+
E_last=IP_val(x, e, x_tilde, m, vol, IB, mu_lame, lam, n, o, contact_area, (x-x_n) /h, mu_lambda, DBC, DBC_target, DBC_stiff[0], h)
33
+
[p, DBC_satisfied] =search_dir(x, e, x_tilde, m, vol, IB, mu_lame, lam, n, o, contact_area, (x-x_n) /h, mu_lambda, is_DBC, DBC, DBC_target, DBC_stiff[0], tol, h)
35
34
while (LA.norm(p, inf) /h>tol) | (sum(DBC_satisfied) !=len(DBC)): # also check whether all DBCs are satisfied
36
35
print('Iteration', iter, ':')
37
36
print('residual =', LA.norm(p, inf) /h)
38
37
39
38
if (LA.norm(p, inf) /h<=tol) & (sum(DBC_satisfied) !=len(DBC)):
40
39
# increase DBC stiffness and recompute energy value record
41
-
DBC_stiff_[0] *=2
42
-
DBC_stiff=DBC_stiff_[0]
43
-
E_last=IP_val(x, e, x_tilde, m, vol, IB, mu_lame, lam, n, o, contact_area, (x-x_n) /h, mu_lambda, DBC, DBC_target, DBC_stiff, h)
40
+
DBC_stiff[0] *=2
41
+
E_last=IP_val(x, e, x_tilde, m, vol, IB, mu_lame, lam, n, o, contact_area, (x-x_n) /h, mu_lambda, DBC, DBC_target, DBC_stiff[0], h)
44
42
45
43
# filter line search
46
44
# ANCHOR: apply_filter
47
45
alpha=min(BarrierEnergy.init_step_size(x, n, o, p), NeoHookeanEnergy.init_step_size(x, e, p)) # avoid interpenetration, tunneling, and inversion
48
46
# ANCHOR_END: apply_filter
49
-
whileIP_val(x+alpha*p, e, x_tilde, m, vol, IB, mu_lame, lam, n, o, contact_area, (x+alpha*p-x_n) /h, mu_lambda, DBC, DBC_target, DBC_stiff, h) >E_last:
47
+
whileIP_val(x+alpha*p, e, x_tilde, m, vol, IB, mu_lame, lam, n, o, contact_area, (x+alpha*p-x_n) /h, mu_lambda, DBC, DBC_target, DBC_stiff[0], h) >E_last:
50
48
alpha/=2
51
49
print('step size =', alpha)
52
50
53
51
x+=alpha*p
54
-
E_last=IP_val(x, e, x_tilde, m, vol, IB, mu_lame, lam, n, o, contact_area, (x-x_n) /h, mu_lambda, DBC, DBC_target, DBC_stiff, h)
55
-
[p, DBC_satisfied] =search_dir(x, e, x_tilde, m, vol, IB, mu_lame, lam, n, o, contact_area, (x-x_n) /h, mu_lambda, is_DBC, DBC, DBC_target, DBC_stiff, tol, h)
52
+
E_last=IP_val(x, e, x_tilde, m, vol, IB, mu_lame, lam, n, o, contact_area, (x-x_n) /h, mu_lambda, DBC, DBC_target, DBC_stiff[0], h)
53
+
[p, DBC_satisfied] =search_dir(x, e, x_tilde, m, vol, IB, mu_lame, lam, n, o, contact_area, (x-x_n) /h, mu_lambda, is_DBC, DBC, DBC_target, DBC_stiff[0], tol, h)
Copy file name to clipboardExpand all lines: 7_self_contact/time_integrator.py
+8-9Lines changed: 8 additions & 9 deletions
Original file line number
Diff line number
Diff line change
@@ -13,7 +13,7 @@
13
13
importFrictionEnergy
14
14
importSpringEnergy
15
15
16
-
defstep_forward(x, e, v, m, vol, IB, mu_lame, lam, n, o, bp, be, contact_area, mu, is_DBC, DBC, DBC_v, DBC_limit, h, tol):
16
+
defstep_forward(x, e, v, m, vol, IB, mu_lame, lam, n, o, bp, be, contact_area, mu, is_DBC, DBC, DBC_v, DBC_limit, DBC_stiff, h, tol):
17
17
x_tilde=x+v*h# implicit Euler predictive position
18
18
x_n=copy.deepcopy(x)
19
19
mu_lambda=BarrierEnergy.compute_mu_lambda(x, n, o, bp, be, contact_area, mu) # compute mu * lambda for each node using x^n
@@ -23,30 +23,29 @@ def step_forward(x, e, v, m, vol, IB, mu_lame, lam, n, o, bp, be, contact_area,
23
23
DBC_target.append(x_n[DBC[i]] +h*DBC_v[i])
24
24
else:
25
25
DBC_target.append(x_n[DBC[i]])
26
-
DBC_stiff=1000# initialize stiffness for DBC springs
27
26
28
27
# Newton loop
29
28
iter=0
30
-
E_last=IP_val(x, e, x_tilde, m, vol, IB, mu_lame, lam, n, o, bp, be, 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, vol, IB, mu_lame, lam, n, o, bp, be, contact_area, (x-x_n) /h, mu_lambda, is_DBC, DBC, DBC_target, DBC_stiff, tol, h)
29
+
E_last=IP_val(x, e, x_tilde, m, vol, IB, mu_lame, lam, n, o, bp, be, contact_area, (x-x_n) /h, mu_lambda, DBC, DBC_target, DBC_stiff[0], h)
30
+
[p, DBC_satisfied] =search_dir(x, e, x_tilde, m, vol, IB, mu_lame, lam, n, o, bp, be, contact_area, (x-x_n) /h, mu_lambda, is_DBC, DBC, DBC_target, DBC_stiff[0], tol, h)
32
31
while (LA.norm(p, inf) /h>tol) | (sum(DBC_satisfied) !=len(DBC)): # also check whether all DBCs are satisfied
33
32
print('Iteration', iter, ':')
34
33
print('residual =', LA.norm(p, inf) /h)
35
34
36
35
if (LA.norm(p, inf) /h<=tol) & (sum(DBC_satisfied) !=len(DBC)):
37
36
# increase DBC stiffness and recompute energy value record
38
-
DBC_stiff*=2
39
-
E_last=IP_val(x, e, x_tilde, m, vol, IB, mu_lame, lam, n, o, bp, be, contact_area, (x-x_n) /h, mu_lambda, DBC, DBC_target, DBC_stiff, h)
37
+
DBC_stiff[0]*=2
38
+
E_last=IP_val(x, e, x_tilde, m, vol, IB, mu_lame, lam, n, o, bp, be, contact_area, (x-x_n) /h, mu_lambda, DBC, DBC_target, DBC_stiff[0], h)
40
39
41
40
# filter line search
42
41
alpha=min(BarrierEnergy.init_step_size(x, n, o, bp, be, p), NeoHookeanEnergy.init_step_size(x, e, p)) # avoid interpenetration, tunneling, and inversion
43
-
whileIP_val(x+alpha*p, e, x_tilde, m, vol, IB, mu_lame, lam, n, o, bp, be, contact_area, (x+alpha*p-x_n) /h, mu_lambda, DBC, DBC_target, DBC_stiff, h) >E_last:
42
+
whileIP_val(x+alpha*p, e, x_tilde, m, vol, IB, mu_lame, lam, n, o, bp, be, contact_area, (x+alpha*p-x_n) /h, mu_lambda, DBC, DBC_target, DBC_stiff[0], h) >E_last:
44
43
alpha/=2
45
44
print('step size =', alpha)
46
45
47
46
x+=alpha*p
48
-
E_last=IP_val(x, e, x_tilde, m, vol, IB, mu_lame, lam, n, o, bp, be, 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, vol, IB, mu_lame, lam, n, o, bp, be, contact_area, (x-x_n) /h, mu_lambda, is_DBC, DBC, DBC_target, DBC_stiff, tol, h)
47
+
E_last=IP_val(x, e, x_tilde, m, vol, IB, mu_lame, lam, n, o, bp, be, contact_area, (x-x_n) /h, mu_lambda, DBC, DBC_target, DBC_stiff[0], h)
48
+
[p, DBC_satisfied] =search_dir(x, e, x_tilde, m, vol, IB, mu_lame, lam, n, o, bp, be, contact_area, (x-x_n) /h, mu_lambda, is_DBC, DBC, DBC_target, DBC_stiff[0], tol, h)
Copy file name to clipboardExpand all lines: 8_self_friction/time_integrator.py
+8-9Lines changed: 8 additions & 9 deletions
Original file line number
Diff line number
Diff line change
@@ -13,7 +13,7 @@
13
13
importFrictionEnergy
14
14
importSpringEnergy
15
15
16
-
defstep_forward(x, e, v, m, vol, IB, mu_lame, lam, n, o, bp, be, contact_area, mu, is_DBC, DBC, DBC_v, DBC_limit, h, tol):
16
+
defstep_forward(x, e, v, m, vol, IB, mu_lame, lam, n, o, bp, be, contact_area, mu, is_DBC, DBC, DBC_v, DBC_limit, DBC_stiff, h, tol):
17
17
x_tilde=x+v*h# implicit Euler predictive position
18
18
x_n=copy.deepcopy(x)
19
19
[mu_lambda, mu_lambda_self] =BarrierEnergy.compute_mu_lambda(x, n, o, bp, be, contact_area, mu) # compute mu * lambda for each node using x^n
@@ -23,30 +23,29 @@ def step_forward(x, e, v, m, vol, IB, mu_lame, lam, n, o, bp, be, contact_area,
23
23
DBC_target.append(x_n[DBC[i]] +h*DBC_v[i])
24
24
else:
25
25
DBC_target.append(x_n[DBC[i]])
26
-
DBC_stiff=1000# initialize stiffness for DBC springs
27
26
28
27
# Newton loop
29
28
iter=0
30
-
E_last=IP_val(x, e, x_tilde, m, vol, IB, mu_lame, lam, n, o, bp, be, contact_area, (x-x_n) /h, mu_lambda, mu_lambda_self, DBC, DBC_target, DBC_stiff, h)
31
-
[p, DBC_satisfied] =search_dir(x, e, x_tilde, m, vol, IB, mu_lame, lam, n, o, bp, be, contact_area, (x-x_n) /h, mu_lambda, mu_lambda_self, is_DBC, DBC, DBC_target, DBC_stiff, tol, h)
29
+
E_last=IP_val(x, e, x_tilde, m, vol, IB, mu_lame, lam, n, o, bp, be, contact_area, (x-x_n) /h, mu_lambda, mu_lambda_self, DBC, DBC_target, DBC_stiff[0], h)
30
+
[p, DBC_satisfied] =search_dir(x, e, x_tilde, m, vol, IB, mu_lame, lam, n, o, bp, be, contact_area, (x-x_n) /h, mu_lambda, mu_lambda_self, is_DBC, DBC, DBC_target, DBC_stiff[0], tol, h)
32
31
while (LA.norm(p, inf) /h>tol) | (sum(DBC_satisfied) !=len(DBC)): # also check whether all DBCs are satisfied
33
32
print('Iteration', iter, ':')
34
33
print('residual =', LA.norm(p, inf) /h)
35
34
36
35
if (LA.norm(p, inf) /h<=tol) & (sum(DBC_satisfied) !=len(DBC)):
37
36
# increase DBC stiffness and recompute energy value record
38
-
DBC_stiff*=2
39
-
E_last=IP_val(x, e, x_tilde, m, vol, IB, mu_lame, lam, n, o, bp, be, contact_area, (x-x_n) /h, mu_lambda, mu_lambda_self, DBC, DBC_target, DBC_stiff, h)
37
+
DBC_stiff[0]*=2
38
+
E_last=IP_val(x, e, x_tilde, m, vol, IB, mu_lame, lam, n, o, bp, be, contact_area, (x-x_n) /h, mu_lambda, mu_lambda_self, DBC, DBC_target, DBC_stiff[0], h)
40
39
41
40
# filter line search
42
41
alpha=min(BarrierEnergy.init_step_size(x, n, o, bp, be, p), NeoHookeanEnergy.init_step_size(x, e, p)) # avoid interpenetration, tunneling, and inversion
43
-
whileIP_val(x+alpha*p, e, x_tilde, m, vol, IB, mu_lame, lam, n, o, bp, be, contact_area, (x+alpha*p-x_n) /h, mu_lambda, mu_lambda_self, DBC, DBC_target, DBC_stiff, h) >E_last:
42
+
whileIP_val(x+alpha*p, e, x_tilde, m, vol, IB, mu_lame, lam, n, o, bp, be, contact_area, (x+alpha*p-x_n) /h, mu_lambda, mu_lambda_self, DBC, DBC_target, DBC_stiff[0], h) >E_last:
44
43
alpha/=2
45
44
print('step size =', alpha)
46
45
47
46
x+=alpha*p
48
-
E_last=IP_val(x, e, x_tilde, m, vol, IB, mu_lame, lam, n, o, bp, be, contact_area, (x-x_n) /h, mu_lambda, mu_lambda_self, DBC, DBC_target, DBC_stiff, h)
49
-
[p, DBC_satisfied] =search_dir(x, e, x_tilde, m, vol, IB, mu_lame, lam, n, o, bp, be, contact_area, (x-x_n) /h, mu_lambda, mu_lambda_self, is_DBC, DBC, DBC_target, DBC_stiff, tol, h)
47
+
E_last=IP_val(x, e, x_tilde, m, vol, IB, mu_lame, lam, n, o, bp, be, contact_area, (x-x_n) /h, mu_lambda, mu_lambda_self, DBC, DBC_target, DBC_stiff[0], h)
48
+
[p, DBC_satisfied] =search_dir(x, e, x_tilde, m, vol, IB, mu_lame, lam, n, o, bp, be, contact_area, (x-x_n) /h, mu_lambda, mu_lambda_self, is_DBC, DBC, DBC_target, DBC_stiff[0], tol, h)
0 commit comments