diff --git a/6_inv_free/simulator.py b/6_inv_free/simulator.py index 28e6e46..6fbde87 100755 --- a/6_inv_free/simulator.py +++ b/6_inv_free/simulator.py @@ -14,7 +14,7 @@ E = 1e5 # Young's modulus nu = 0.4 # Poisson's ratio # ANCHOR_END: lame_param -n_seg = 4 # num of segments per side of the square +n_seg = 4 # num of segments per side of the square h = 0.01 # time step size in s DBC = [(n_seg + 1) * (n_seg + 1)] # dirichlet node index DBC_v = [np.array([0.0, -0.5])] # dirichlet node velocity @@ -50,6 +50,7 @@ resolution = np.array([900, 900]) offset = resolution / 2 scale = 200 +DBC_stiff = [1000] def screen_projection(x): return [offset[0] + scale * x[0], resolution[1] - (offset[1] + scale * x[1])] @@ -84,7 +85,7 @@ def screen_projection(x): pygame.display.flip() # flip the display # step forward simulation and wait for screen refresh - [x, v] = time_integrator.step_forward(x, e, v, m, vol, IB, mu_lame, lam, ground_n, ground_o, contact_area, mu, is_DBC, DBC, DBC_v, DBC_limit, h, 1e-2) + [x, v] = time_integrator.step_forward(x, e, v, m, vol, IB, mu_lame, lam, ground_n, ground_o, contact_area, mu, is_DBC, DBC, DBC_v, DBC_limit,DBC_stiff, h, 1e-2) time_step += 1 pygame.time.wait(int(h * 1000)) square_mesh.write_to_file(time_step, x, e) diff --git a/6_inv_free/time_integrator.py b/6_inv_free/time_integrator.py index 7e43a41..4b5fc4e 100644 --- a/6_inv_free/time_integrator.py +++ b/6_inv_free/time_integrator.py @@ -13,7 +13,10 @@ import FrictionEnergy import SpringEnergy -def step_forward(x, e, v, m, vol, IB, mu_lame, lam, n, o, contact_area, mu, is_DBC, DBC, DBC_v, DBC_limit, h, tol): + + + +def step_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): x_tilde = x + v * h # implicit Euler predictive position x_n = copy.deepcopy(x) mu_lambda = BarrierEnergy.compute_mu_lambda(x, n, o, contact_area, mu) # compute mu * lambda for each node using x^n @@ -23,7 +26,7 @@ def step_forward(x, e, v, m, vol, IB, mu_lame, lam, n, o, contact_area, mu, is_D DBC_target.append(x_n[DBC[i]] + h * DBC_v[i]) else: DBC_target.append(x_n[DBC[i]]) - DBC_stiff = 1000 # initialize stiffness for DBC springs + DBC_stiff = DBC_stiff_[0] # Newton loop iter = 0 @@ -35,7 +38,8 @@ def step_forward(x, e, v, m, vol, IB, mu_lame, lam, n, o, contact_area, mu, is_D if (LA.norm(p, inf) / h <= tol) & (sum(DBC_satisfied) != len(DBC)): # increase DBC stiffness and recompute energy value record - DBC_stiff *= 2 + DBC_stiff_[0] *= 2 + DBC_stiff = DBC_stiff_[0] 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) # filter line search