@@ -14,10 +14,10 @@ template <typename T, int dim>
1414struct MovDirichletSimulator <T, dim>::Impl
1515{
1616 int n_seg;
17- T h, rho, side_len, initial_stretch, m, tol, mu;
17+ T h, rho, side_len, initial_stretch, m, tol, mu, DBC_stiff ;
1818 int resolution = 900 , scale = 200 , offset = resolution / 2 , radius = 5 ;
19- std::vector<T> x, x_tilde, v, k, l2;
20- std::vector<int > e;
19+ std::vector<T> x, x_tilde, v, k, l2, DBC_limit, DBC_v, DBC_target ;
20+ std::vector<int > e, DBC, DBC_satisfied ;
2121 DeviceBuffer<int > device_DBC;
2222 DeviceBuffer<T> device_contact_area;
2323 sf::RenderWindow window;
@@ -32,6 +32,7 @@ struct MovDirichletSimulator<T, dim>::Impl
3232 void update_x_tilde (const DeviceBuffer<T> &new_x_tilde);
3333 void update_v (const DeviceBuffer<T> &new_v);
3434 void update_DBC_target ();
35+ void update_DBC_stiff (T new_DBC_stiff);
3536 T IP_val ();
3637 void step_forward ();
3738 void draw ();
@@ -61,10 +62,18 @@ template <typename T, int dim>
6162MovDirichletSimulator<T, dim>::Impl::Impl(T rho, T side_len, T initial_stretch, T K, T h_, T tol_, T mu_, int n_seg) : tol(tol_), h(h_), mu(mu_), window(sf::VideoMode(resolution, resolution), " MovDirichletSimulator" )
6263{
6364 generate (side_len, n_seg, x, e);
64- std::vector<int > DBC (x.size () / dim, 0 );
65+ DBC.push_back ((n_seg + 1 ) * (n_seg + 1 ));
66+ DBC_target.resize (DBC.size () * dim);
67+ DBC_limit.push_back (0 );
68+ DBC_limit.push_back (-0.6 );
69+ DBC_v.push_back (0 );
70+ DBC_v.push_back (-0.5 );
71+ DBC_stiff = 10 ;
72+ x.push_back (0 );
73+ x.push_back (side_len * 0.6 );
6574 std::vector<T> contact_area (x.size () / dim, side_len / n_seg);
6675 std::vector<T> ground_n (dim);
67- ground_n[0 ] = 0.1 , ground_n[1 ] = 1 ;
76+ ground_n[0 ] = 0 , ground_n[1 ] = 1 ;
6877 T n_norm = ground_n[0 ] * ground_n[0 ] + ground_n[1 ] * ground_n[1 ];
6978 n_norm = sqrt (n_norm);
7079 for (int i = 0 ; i < dim; i++)
@@ -94,7 +103,7 @@ MovDirichletSimulator<T, dim>::Impl::Impl(T rho, T side_len, T initial_stretch,
94103 gravityenergy = GravityEnergy<T, dim>(N, m);
95104 barrierenergy = BarrierEnergy<T, dim>(x, ground_n, ground_o, contact_area);
96105 frictionenergy = FrictionEnergy<T, dim>(v, h, ground_n);
97- springenergy = SpringEnergy<T, dim>(x, std::vector<T>(N, m), DBC, std::vector<T>(N * dim, 0 ), std::vector<T>(N * dim, 0 ), 0 , h );
106+ springenergy = SpringEnergy<T, dim>(x, std::vector<T>(N, m), DBC, DBC_stiff );
98107 DeviceBuffer<T> x_device (x);
99108 update_x (x_device);
100109 device_DBC = DeviceBuffer<int >(DBC);
@@ -132,15 +141,22 @@ void MovDirichletSimulator<T, dim>::Impl::step_forward()
132141 DeviceBuffer<T> x_tilde (x.size ()); // Predictive position
133142 update_x_tilde (add_vector<T>(x, v, 1 , h));
134143 frictionenergy.update_mu_lambda (barrierenergy.compute_mu_lambda (mu));
144+ update_DBC_target ();
145+ update_DBC_stiff (10 );
135146 DeviceBuffer<T> x_n = x; // Copy current positions to x_n
136147 update_v (add_vector<T>(x, x_n, 1 / h, -1 / h));
137148 int iter = 0 ;
138149 T E_last = IP_val ();
139150 DeviceBuffer<T> p = search_direction ();
140151 T residual = max_vector (p) / h;
141152 // std::cout << "Initial residual " << residual << "\n";
142- while (residual > tol)
153+ while (residual > tol || DBC_satisfied. back () != 1 ) // use last one for simplisity, should check all
143154 {
155+ if (residual <= tol && DBC_satisfied.back () != 1 )
156+ {
157+ update_DBC_stiff (DBC_stiff * 2 );
158+ E_last = IP_val ();
159+ }
144160 // Line search
145161 T alpha = barrierenergy.init_step_size (p);
146162 DeviceBuffer<T> x0 = x;
@@ -195,7 +211,35 @@ void MovDirichletSimulator<T, dim>::Impl::update_v(const DeviceBuffer<T> &new_v)
195211template <typename T, int dim>
196212void MovDirichletSimulator<T, dim>::Impl::update_DBC_target()
197213{
198- springenergy.update_DBC_target ();
214+ for (int i = 0 ; i < DBC.size (); i++)
215+ {
216+ T diff = 0 ;
217+ for (int d = 0 ; d < dim; d++)
218+ {
219+ diff += (DBC_limit[i * dim + d] - x[DBC[i] * dim + d]) * DBC_v[i * dim + d];
220+ }
221+ if (diff > 0 )
222+ {
223+ for (int d = 0 ; d < dim; d++)
224+ {
225+ DBC_target[i * dim + d] = x[DBC[i] * dim + d] + h * DBC_v[i * dim + d];
226+ }
227+ }
228+ else
229+ {
230+ for (int d = 0 ; d < dim; d++)
231+ {
232+ DBC_target[i * dim + d] = x[DBC[i] * dim + d];
233+ }
234+ }
235+ }
236+ springenergy.update_DBC_target (DBC_target);
237+ }
238+ template <typename T, int dim>
239+ void MovDirichletSimulator<T, dim>::Impl::update_DBC_stiff(T new_DBC_stiff)
240+ {
241+ DBC_stiff = new_DBC_stiff;
242+ springenergy.update_k (new_DBC_stiff);
199243}
200244template <typename T, int dim>
201245void MovDirichletSimulator<T, dim>::Impl::draw()
@@ -212,28 +256,27 @@ void MovDirichletSimulator<T, dim>::Impl::draw()
212256 }
213257
214258 // Draw masses as circles
215- for (int i = 0 ; i < x.size () / dim; ++i)
259+ for (int i = 0 ; i < ( x.size () - 1 ) / dim; ++i)
216260 {
217261 sf::CircleShape circle (radius); // Set a fixed radius for each mass
218262 circle.setFillColor (sf::Color::Red);
219263 circle.setPosition (screen_projection_x (x[i * dim]) - radius, screen_projection_y (x[i * dim + 1 ]) - radius); // Center the circle on the mass
220264 window.draw (circle);
221265 }
222-
223266 window.display (); // Display the rendered frame
224267}
225268
226269template <typename T, int dim>
227270T MovDirichletSimulator<T, dim>::Impl::IP_val()
228271{
229272
230- return inertialenergy.val () + (massspringenergy.val () + gravityenergy.val () + barrierenergy.val () + frictionenergy.val ()) * h * h;
273+ return inertialenergy.val () + (massspringenergy.val () + gravityenergy.val () + barrierenergy.val () + frictionenergy.val ()) * h * h + springenergy. val () ;
231274}
232275
233276template <typename T, int dim>
234277DeviceBuffer<T> MovDirichletSimulator<T, dim>::Impl::IP_grad()
235278{
236- return add_vector<T>(add_vector<T>(add_vector<T>(add_vector<T>(inertialenergy.grad (), massspringenergy.grad (), 1.0 , h * h), gravityenergy.grad (), 1.0 , h * h), barrierenergy.grad (), 1.0 , h * h), frictionenergy.grad (), 1.0 , h * h);
279+ return add_vector<T>(add_vector<T>(add_vector<T>(add_vector<T>(add_vector<T>( inertialenergy.grad (), massspringenergy.grad (), 1.0 , h * h), gravityenergy.grad (), 1.0 , h * h), barrierenergy.grad (), 1.0 , h * h), frictionenergy.grad (), 1.0 , h * h), springenergy. grad (), 1.0 , 1.0 );
237280}
238281
239282template <typename T, int dim>
@@ -246,6 +289,8 @@ DeviceTripletMatrix<T, 1> MovDirichletSimulator<T, dim>::Impl::IP_hess()
246289 hess = add_triplet<T>(hess, barrier_hess, 1.0 , h * h);
247290 DeviceTripletMatrix<T, 1 > friction_hess = frictionenergy.hess ();
248291 hess = add_triplet<T>(hess, friction_hess, 1.0 , h * h);
292+ DeviceTripletMatrix<T, 1 > spring_hess = springenergy.hess ();
293+ hess = add_triplet<T>(hess, spring_hess, 1.0 , 1.0 );
249294 return hess;
250295}
251296template <typename T, int dim>
@@ -255,7 +300,21 @@ DeviceBuffer<T> MovDirichletSimulator<T, dim>::Impl::search_direction()
255300 dir.resize (x.size ());
256301 DeviceBuffer<T> grad = IP_grad ();
257302 DeviceTripletMatrix<T, 1 > hess = IP_hess ();
258- search_dir<T, dim>(grad, hess, dir, device_DBC);
303+ // check whether each DBC is satisfied
304+ DBC_satisfied.resize (x.size () / dim, 0 );
305+ for (int i = 0 ; i < DBC.size (); i++)
306+ {
307+ T diff = 0 ;
308+ for (int d = 0 ; d < dim; d++)
309+ {
310+ diff += (x[DBC[i] * dim + d] - DBC_target[i * dim + d]) * (x[DBC[i] * dim + d] - DBC_target[i * dim + d]);
311+ }
312+ if (diff / h < tol)
313+ {
314+ DBC_satisfied[DBC[i]] = 1 ;
315+ }
316+ }
317+ search_dir<T, dim>(grad, hess, dir, device_DBC, DBC_target, DBC_satisfied);
259318 return dir;
260319}
261320
0 commit comments