Skip to content

Commit e1ccc79

Browse files
RoushelfyRoushelfy
committed
4_friction
1 parent f32f249 commit e1ccc79

8 files changed

Lines changed: 68 additions & 67 deletions

File tree

‎simulators/4_friction/include/BarrierEnergy.h‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -19,7 +19,7 @@ class BarrierEnergy
1919
T val(); // Calculate the value of the energy
2020
const DeviceBuffer<T> &grad(); // Calculate the gradient of the energy
2121
const DeviceTripletMatrix<T, 1> &hess(); // Calculate the Hessian matrix of the energy
22-
const DeviceBuffer<T> &compute_mu_lambda();
22+
const DeviceBuffer<T> compute_mu_lambda(T mu);
2323
T init_step_size(const DeviceBuffer<T> &p); // Calculate the initial step size for the line search
2424

2525
private:

‎simulators/4_friction/include/FrictionEnergy.h‎

Lines changed: 5 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -9,22 +9,23 @@ template <typename T, int dim>
99
class FrictionEnergy
1010
{
1111
public:
12-
FrictionEnergy(const std::vector<T> &v, const std::vector<T> &mu_lambda, T hhat, const Eigen::Matrix<T, dim, 1> &n);
12+
FrictionEnergy(const std::vector<T> &v, T hhat, const std::vector<T> &n);
1313
FrictionEnergy();
1414
~FrictionEnergy();
1515
FrictionEnergy(FrictionEnergy &&rhs);
1616
FrictionEnergy(const FrictionEnergy &rhs);
1717
FrictionEnergy &operator=(FrictionEnergy &&rhs);
1818

1919
void update_v(const DeviceBuffer<T> &v);
20+
void update_mu_lambda(const DeviceBuffer<T> &mu_lambda);
2021
T val(); // Calculate the value of the energy
2122
const DeviceBuffer<T> &grad(); // Calculate the gradient of the energy
2223
const DeviceTripletMatrix<T, 1> &hess(); // Calculate the Hessian matrix of the energy
2324

2425
private:
2526
struct Impl;
2627
std::unique_ptr<Impl> pimpl_;
27-
T f0(T vbarnorm, T Epsv, T hhat);
28-
T f1_div_vbarnorm(T vbarnorm, T Epsv);
29-
T f_hess_term(T vbarnorm, T Epsv);
28+
T __device__ f0(T vbarnorm, T Epsv, T hhat);
29+
T __device__ f1_div_vbarnorm(T vbarnorm, T Epsv);
30+
T __device__ f_hess_term(T vbarnorm, T Epsv);
3031
};

‎simulators/4_friction/include/simulator.h‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -12,7 +12,7 @@ class FrictionSimulator
1212
~FrictionSimulator();
1313
FrictionSimulator(FrictionSimulator &&rhs);
1414
FrictionSimulator &operator=(FrictionSimulator &&rhs);
15-
FrictionSimulator(T rho, T side_len, T initial_stretch, T K, T h, T tol, int n_seg);
15+
FrictionSimulator(T rho, T side_len, T initial_stretch, T K, T h, T tol, T mu, int n_seg);
1616
void run();
1717

1818
private:

‎simulators/4_friction/src/BarrierEnergy.cu‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -184,13 +184,13 @@ template class BarrierEnergy<double, 2>;
184184
template class BarrierEnergy<double, 3>;
185185

186186
template <typename T, int dim>
187-
const DeviceBuffer<T> &BarrierEnergy<T, dim>::compute_mu_lambda(T mu)
187+
const DeviceBuffer<T> BarrierEnergy<T, dim>::compute_mu_lambda(T mu)
188188
{
189189
auto &device_x = pimpl_->device_x;
190190
auto &device_n = pimpl_->device_n;
191191
auto &device_o = pimpl_->device_o;
192192
auto &device_contact_area = pimpl_->device_contact_area;
193-
int N = device_x.size();
193+
int N = device_x.size() / dim;
194194
DeviceBuffer<T> device_mu_lambda(N);
195195
device_mu_lambda.fill(0);
196196
ParallelFor(256)

‎simulators/4_friction/src/FrictionEnergy.cu‎

Lines changed: 37 additions & 51 deletions
Original file line numberDiff line numberDiff line change
@@ -36,14 +36,14 @@ FrictionEnergy<T, dim>::FrictionEnergy(const FrictionEnergy<T, dim> &rhs)
3636
: pimpl_{std::make_unique<Impl>(*rhs.pimpl_)} {}
3737

3838
template <typename T, int dim>
39-
FrictionEnergy<T, dim>::FrictionEnergy(const std::vector<T> &v, const std::vector<T> &mu_lambda, T hhat, const Eigen::Matrix<T, dim, 1> &n)
39+
FrictionEnergy<T, dim>::FrictionEnergy(const std::vector<T> &v, T hhat, const std::vector<T> &n)
4040
: pimpl_{std::make_unique<Impl>()}
4141
{
4242
pimpl_->N = v.size() / dim;
4343
pimpl_->device_v.copy_from(v);
44-
pimpl_->device_mu_lambda.copy_from(mu_lambda);
44+
pimpl_->device_mu_lambda.resize(pimpl_->N);
4545
pimpl_->hhat = hhat;
46-
pimpl_->n = n;
46+
pimpl_->n = Eigen::Map<const Eigen::Matrix<T, dim, 1>>(n.data());
4747
pimpl_->device_grad.resize(pimpl_->N * dim);
4848
pimpl_->device_hess.resize_triplets(pimpl_->N * dim * dim);
4949
pimpl_->device_hess.reshape(v.size(), v.size());
@@ -54,9 +54,14 @@ void FrictionEnergy<T, dim>::update_v(const DeviceBuffer<T> &v)
5454
{
5555
pimpl_->device_v.view().copy_from(v);
5656
}
57+
template <typename T, int dim>
58+
void FrictionEnergy<T, dim>::update_mu_lambda(const DeviceBuffer<T> &mu_lambda)
59+
{
60+
pimpl_->device_mu_lambda.view().copy_from(mu_lambda);
61+
}
5762

5863
template <typename T, int dim>
59-
T FrictionEnergy<T, dim>::f0(T vbarnorm, T Epsv, T hhat)
64+
T __device__ FrictionEnergy<T, dim>::f0(T vbarnorm, T Epsv, T hhat)
6065
{
6166
if (vbarnorm >= Epsv)
6267
{
@@ -65,13 +70,13 @@ T FrictionEnergy<T, dim>::f0(T vbarnorm, T Epsv, T hhat)
6570
else
6671
{
6772
T vbarnormhhat = vbarnorm * hhat;
68-
T Epsvhhat = Epsv * hhat;
73+
T epsvhhat = Epsv * hhat;
6974
return vbarnormhhat * vbarnormhhat * (-vbarnormhhat / 3.0 + epsvhhat) / (epsvhhat * epsvhhat) + epsvhhat / 3.0;
7075
}
7176
}
7277

7378
template <typename T, int dim>
74-
T FrictionEnergy<T, dim>::f1_div_vbarnorm(T vbarnorm, T Epsv)
79+
T __device__ FrictionEnergy<T, dim>::f1_div_vbarnorm(T vbarnorm, T Epsv)
7580
{
7681
if (vbarnorm >= Epsv)
7782
{
@@ -84,7 +89,7 @@ T FrictionEnergy<T, dim>::f1_div_vbarnorm(T vbarnorm, T Epsv)
8489
}
8590

8691
template <typename T, int dim>
87-
T FrictionEnergy<T, dim>::f_hess_term(T vbarnorm, T Epsv)
92+
T __device__ FrictionEnergy<T, dim>::f_hess_term(T vbarnorm, T Epsv)
8893
{
8994
if (vbarnorm >= Epsv)
9095
{
@@ -97,7 +102,7 @@ T FrictionEnergy<T, dim>::f_hess_term(T vbarnorm, T Epsv)
97102
}
98103

99104
template <typename T, int dim>
100-
T Energy<T, dim>::val()
105+
T FrictionEnergy<T, dim>::val()
101106
{
102107
auto &device_v = pimpl_->device_v;
103108
auto &device_mu_lambda = pimpl_->device_mu_lambda;
@@ -111,7 +116,11 @@ T Energy<T, dim>::val()
111116
Eigen::Matrix<T, dim, dim> T_mat = Eigen::Matrix<T, dim, dim>::Identity() - n * n.transpose();
112117
if (device_mu_lambda(i) > 0)
113118
{
114-
Eigen::Matrix<T, dim, 1> v = device_v.segment(i * dim, dim);
119+
Eigen::Matrix<T, dim, 1> v;
120+
for (int j = 0; j < dim; ++j)
121+
{
122+
v(j) = device_v(i * dim + j);
123+
}
115124
Eigen::Matrix<T, dim, 1> vbar = T_mat * v;
116125
T vbarnorm = vbar.norm();
117126
T val = f0(vbarnorm, epsv, hhat);
@@ -123,7 +132,7 @@ T Energy<T, dim>::val()
123132
}
124133

125134
template <typename T, int dim>
126-
const DeviceBuffer<T> &Energy<T, dim>::grad()
135+
const DeviceBuffer<T> &FrictionEnergy<T, dim>::grad()
127136
{
128137
auto &device_v = pimpl_->device_v;
129138
auto &device_mu_lambda = pimpl_->device_mu_lambda;
@@ -138,7 +147,11 @@ const DeviceBuffer<T> &Energy<T, dim>::grad()
138147
Eigen::Matrix<T, dim, dim> T_mat = Eigen::Matrix<T, dim, dim>::Identity() - n * n.transpose();
139148
if (device_mu_lambda(i) > 0)
140149
{
141-
Eigen::Matrix<T, dim, 1> v = device_v.segment(i * dim, dim);
150+
Eigen::Matrix<T, dim, 1> v;
151+
for (int j = 0; j < dim; ++j)
152+
{
153+
v(j) = device_v(i * dim + j);
154+
}
142155
Eigen::Matrix<T, dim, 1> vbar = T_mat * v;
143156
T vbarnorm = vbar.norm();
144157
T grad_factor = f1_div_vbarnorm(vbarnorm, epsv);
@@ -153,7 +166,6 @@ const DeviceBuffer<T> &Energy<T, dim>::grad()
153166

154167
return device_grad;
155168
}
156-
157169
template <typename T, int dim>
158170
const DeviceTripletMatrix<T, 1> &FrictionEnergy<T, dim>::hess()
159171
{
@@ -168,64 +180,38 @@ const DeviceTripletMatrix<T, 1> &FrictionEnergy<T, dim>::hess()
168180
int N = device_v.size() / dim;
169181
ParallelFor(256).apply(N, [device_v = device_v.cviewer(), device_mu_lambda = device_mu_lambda.cviewer(), device_hess_val = device_hess_val.viewer(), device_hess_row_idx = device_hess_row_idx.viewer(), device_hess_col_idx = device_hess_col_idx.viewer(), hhat, n, N, this] __device__(int i) mutable
170182
{
171-
T T_mat[dim][dim];
172-
for(int r = 0; r < dim; ++r)
173-
{
174-
for(int c = 0; c < dim; ++c)
175-
{
176-
T_mat[r][c] = (r == c ? 1.0 : 0.0) - n(r) * n(c);
177-
}
178-
}
183+
Eigen::Matrix<T, dim, dim> T_mat = Eigen::Matrix<T, dim, dim>::Identity() - n * n.transpose();
179184
if (device_mu_lambda(i) > 0)
180185
{
181-
T vbar[dim] = {0};
186+
Eigen::Matrix<T, dim, 1> v;
182187
for (int j = 0; j < dim; ++j)
183188
{
184-
for (int k = 0; k < dim; ++k)
185-
{
186-
vbar[j] += T_mat[j][k] * device_v(i * dim + k);
187-
}
188-
}
189-
T vbarnorm = 0;
190-
for (int j = 0; j < dim; ++j)
191-
{
192-
vbarnorm += vbar[j] * vbar[j];
193-
}
194-
vbarnorm = sqrt(vbarnorm);
195-
T inner_term[dim][dim];
196-
for(int r = 0; r < dim; ++r)
197-
{
198-
for(int c = 0; c < dim; ++c)
199-
{
200-
inner_term[r][c] = (r == c ? 1.0 : 0.0);
201-
}
189+
v(j) = device_v(i * dim + j);
202190
}
203-
T hess_factor = f_hess_term(vbarnorm, epsv);
191+
Eigen::Matrix<T, dim, 1> vbar = T_mat * v;
192+
T vbarnorm = vbar.norm();
193+
Eigen::Matrix<T, dim, dim> inner_term = Eigen::Matrix<T, dim, dim>::Identity() * f1_div_vbarnorm(vbarnorm, epsv);
204194
if (vbarnorm != 0)
205195
{
206-
for(int r = 0; r < dim; ++r)
207-
{
208-
for(int c = 0; c < dim; ++c)
209-
{
210-
inner_term[r][c] += hess_factor / vbarnorm * vbar[r] * vbar[c];
211-
}
212-
}
196+
inner_term += f_hess_term(vbarnorm, epsv) / vbarnorm * vbar * vbar.transpose();
213197
}
214-
for(int j = 0; j < dim; ++j)
198+
Eigen::Matrix<T, dim, dim> local_hess;
199+
make_PSD(inner_term, local_hess);
200+
local_hess = device_mu_lambda(i) * T_mat * local_hess * T_mat.transpose() / hhat;
201+
for (int j = 0; j < dim; ++j)
215202
{
216-
for(int k = 0; k < dim; ++k)
203+
for (int k = 0; k < dim; ++k)
217204
{
218205
int idx = i * dim * dim + j * dim + k;
219206
device_hess_row_idx(idx) = i * dim + j;
220207
device_hess_col_idx(idx) = i * dim + k;
221-
device_hess_val(idx) = device_mu_lambda(i) * inner_term[j][k] / hhat;
208+
device_hess_val(idx) = local_hess(j, k);
222209
}
223210
}
224211
} })
225212
.wait();
226213
return device_hess;
227214
}
228-
229215
template class FrictionEnergy<float, 2>;
230216
template class FrictionEnergy<float, 3>;
231217
template class FrictionEnergy<double, 2>;

‎simulators/4_friction/src/device_uti.cu‎

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -38,6 +38,10 @@ void __device__ make_PSD(const Eigen::Matrix<T, Size, Size> &hess, Eigen::Matrix
3838
PSD = V * lamDiag * VT;
3939
}
4040

41+
template void __device__ make_PSD<float, 2>(const Eigen::Matrix<float, 2, 2> &hess, Eigen::Matrix<float, 2, 2> &PSD);
42+
template void __device__ make_PSD<double, 2>(const Eigen::Matrix<double, 2, 2> &hess, Eigen::Matrix<double, 2, 2> &PSD);
43+
template void __device__ make_PSD<float, 3>(const Eigen::Matrix<float, 3, 3> &hess, Eigen::Matrix<float, 3, 3> &PSD);
44+
template void __device__ make_PSD<double, 3>(const Eigen::Matrix<double, 3, 3> &hess, Eigen::Matrix<double, 3, 3> &PSD);
4145
template void __device__ make_PSD<float, 4>(const Eigen::Matrix<float, 4, 4> &hess, Eigen::Matrix<float, 4, 4> &PSD);
4246
template void __device__ make_PSD<double, 4>(const Eigen::Matrix<double, 4, 4> &hess, Eigen::Matrix<double, 4, 4> &PSD);
4347
template void __device__ make_PSD<float, 6>(const Eigen::Matrix<float, 6, 6> &hess, Eigen::Matrix<float, 6, 6> &PSD);

‎simulators/4_friction/src/main.cpp‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -2,8 +2,8 @@
22

33
int main()
44
{
5-
float rho = 1000, k = 4e4, initial_stretch = 1, n_seg = 10, h = 0.01, side_len = 1, tol = 0.01;
5+
float rho = 1000, k = 4e4, initial_stretch = 1, n_seg = 8, h = 0.01, side_len = 1, tol = 0.01, mu = 0.11;
66
// printf("Running mass-spring simulator with parameters: rho = %f, k = %f, initial_stretch = %f, n_seg = %d, h = %f, side_len = %f, tol = %f\n", rho, k, initial_stretch, n_seg, h, side_len, tol);
7-
FrictionSimulator<float, 2> simulator(rho, side_len, initial_stretch, k, h, tol, n_seg);
7+
FrictionSimulator<float, 2> simulator(rho, side_len, initial_stretch, k, h, tol, mu, n_seg);
88
simulator.run();
99
}

‎simulators/4_friction/src/simulator.cu‎

Lines changed: 16 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -4,6 +4,7 @@
44
#include "MassSpringEnergy.h"
55
#include "GravityEnergy.h"
66
#include "BarrierEnergy.h"
7+
#include "FrictionEnergy.h"
78
#include <muda/muda.h>
89
#include <muda/container.h>
910
#include "uti.h"
@@ -12,7 +13,7 @@ template <typename T, int dim>
1213
struct FrictionSimulator<T, dim>::Impl
1314
{
1415
int n_seg;
15-
T h, rho, side_len, initial_stretch, m, tol;
16+
T h, rho, side_len, initial_stretch, m, tol, mu;
1617
int resolution = 900, scale = 200, offset = resolution / 2, radius = 5;
1718
std::vector<T> x, x_tilde, v, k, l2;
1819
std::vector<int> e;
@@ -23,7 +24,8 @@ struct FrictionSimulator<T, dim>::Impl
2324
MassSpringEnergy<T, dim> massspringenergy;
2425
GravityEnergy<T, dim> gravityenergy;
2526
BarrierEnergy<T, dim> barrierenergy;
26-
Impl(T rho, T side_len, T initial_stretch, T K, T h_, T tol_, int n_seg);
27+
FrictionEnergy<T, dim> frictionenergy;
28+
Impl(T rho, T side_len, T initial_stretch, T K, T h_, T tol_, T mu_, int n_seg);
2729
void update_x(const DeviceBuffer<T> &new_x);
2830
void update_x_tilde(const DeviceBuffer<T> &new_x_tilde);
2931
void update_v(const DeviceBuffer<T> &new_v);
@@ -49,11 +51,11 @@ template <typename T, int dim>
4951
FrictionSimulator<T, dim> &FrictionSimulator<T, dim>::operator=(FrictionSimulator<T, dim> &&rhs) = default;
5052

5153
template <typename T, int dim>
52-
FrictionSimulator<T, dim>::FrictionSimulator(T rho, T side_len, T initial_stretch, T K, T h_, T tol_, int n_seg) : pimpl_{std::make_unique<Impl>(rho, side_len, initial_stretch, K, h_, tol_, n_seg)}
54+
FrictionSimulator<T, dim>::FrictionSimulator(T rho, T side_len, T initial_stretch, T K, T h_, T tol_, T mu_, int n_seg) : pimpl_{std::make_unique<Impl>(rho, side_len, initial_stretch, K, h_, tol_, mu_, n_seg)}
5355
{
5456
}
5557
template <typename T, int dim>
56-
FrictionSimulator<T, dim>::Impl::Impl(T rho, T side_len, T initial_stretch, T K, T h_, T tol_, int n_seg) : tol(tol_), h(h_), window(sf::VideoMode(resolution, resolution), "FrictionSimulator")
58+
FrictionSimulator<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), "FrictionSimulator")
5759
{
5860
generate(side_len, n_seg, x, e);
5961
std::vector<int> DBC(x.size() / dim, 0);
@@ -88,6 +90,7 @@ FrictionSimulator<T, dim>::Impl::Impl(T rho, T side_len, T initial_stretch, T K,
8890
massspringenergy = MassSpringEnergy<T, dim>(x, e, l2, k);
8991
gravityenergy = GravityEnergy<T, dim>(N, m);
9092
barrierenergy = BarrierEnergy<T, dim>(x, ground_n, ground_o, contact_area);
93+
frictionenergy = FrictionEnergy<T, dim>(v, h, ground_n);
9194
DeviceBuffer<T> x_device(x);
9295
update_x(x_device);
9396
device_DBC = DeviceBuffer<int>(DBC);
@@ -124,7 +127,9 @@ void FrictionSimulator<T, dim>::Impl::step_forward()
124127
{
125128
DeviceBuffer<T> x_tilde(x.size()); // Predictive position
126129
update_x_tilde(add_vector<T>(x, v, 1, h));
130+
frictionenergy.update_mu_lambda(barrierenergy.compute_mu_lambda(mu));
127131
DeviceBuffer<T> x_n = x; // Copy current positions to x_n
132+
update_v(add_vector<T>(x, x_n, 1 / h, -1 / h));
128133
int iter = 0;
129134
T E_last = IP_val();
130135
DeviceBuffer<T> p = search_direction();
@@ -136,10 +141,12 @@ void FrictionSimulator<T, dim>::Impl::step_forward()
136141
T alpha = barrierenergy.init_step_size(p);
137142
DeviceBuffer<T> x0 = x;
138143
update_x(add_vector<T>(x0, p, 1.0, alpha));
144+
update_v(add_vector<T>(x, x_n, 1 / h, -1 / h));
139145
while (IP_val() > E_last)
140146
{
141147
alpha /= 2;
142148
update_x(add_vector<T>(x0, p, 1.0, alpha));
149+
update_v(add_vector<T>(x, x_n, 1 / h, -1 / h));
143150
}
144151
std::cout << "step size = " << alpha << "\n";
145152
E_last = IP_val();
@@ -178,6 +185,7 @@ void FrictionSimulator<T, dim>::Impl::update_x_tilde(const DeviceBuffer<T> &new_
178185
template <typename T, int dim>
179186
void FrictionSimulator<T, dim>::Impl::update_v(const DeviceBuffer<T> &new_v)
180187
{
188+
frictionenergy.update_v(new_v);
181189
new_v.copy_to(v);
182190
}
183191
template <typename T, int dim>
@@ -210,13 +218,13 @@ template <typename T, int dim>
210218
T FrictionSimulator<T, dim>::Impl::IP_val()
211219
{
212220

213-
return inertialenergy.val() + (massspringenergy.val() + gravityenergy.val() + barrierenergy.val()) * h * h;
221+
return inertialenergy.val() + (massspringenergy.val() + gravityenergy.val() + barrierenergy.val() + frictionenergy.val()) * h * h;
214222
}
215223

216224
template <typename T, int dim>
217225
DeviceBuffer<T> FrictionSimulator<T, dim>::Impl::IP_grad()
218226
{
219-
return 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);
227+
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);
220228
}
221229

222230
template <typename T, int dim>
@@ -227,6 +235,8 @@ DeviceTripletMatrix<T, 1> FrictionSimulator<T, dim>::Impl::IP_hess()
227235
DeviceTripletMatrix<T, 1> hess = add_triplet<T>(inertial_hess, massspring_hess, 1.0, h * h);
228236
DeviceTripletMatrix<T, 1> barrier_hess = barrierenergy.hess();
229237
hess = add_triplet<T>(hess, barrier_hess, 1.0, h * h);
238+
DeviceTripletMatrix<T, 1> friction_hess = frictionenergy.hess();
239+
hess = add_triplet<T>(hess, friction_hess, 1.0, h * h);
230240
return hess;
231241
}
232242
template <typename T, int dim>

0 commit comments

Comments
 (0)