Skip to content

Commit 3ca4736

Browse files
RoushelfyRoushelfy
committed
7
1 parent cdc96fb commit 3ca4736

8 files changed

Lines changed: 78 additions & 48 deletions

File tree

‎simulators/6_inv_free/src/BarrierEnergy.cu‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -241,7 +241,7 @@ T BarrierEnergy<T, dim>::init_step_size(const DeviceBuffer<T> &p)
241241
T p_n = 0;
242242
for (int j = 0; j < dim; j++)
243243
{
244-
p_n += p(i * dim + j) * device_n_ceil(j);
244+
p_n += (p(i * dim + j)-p((N-1)*dim+j)) * device_n_ceil(j);
245245
}
246246
if (p_n < 0)
247247
{

‎simulators/6_inv_free/src/NeoHookeanEnergy.cu‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -96,7 +96,8 @@ T NeoHookeanEnergy<T, dim>::val()
9696
NeoHookeanEnergyVal(E, Mu, Lambda,X,device_IB(i),device_vol(i));
9797
device_val(i) = E; })
9898
.wait();
99-
return devicesum(device_val);
99+
T val = devicesum(device_val);
100+
return val;
100101
}
101102

102103
template <typename T, int dim>

‎simulators/6_inv_free/src/main.cpp‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -6,7 +6,7 @@ int main()
66
double nu = 0.4, E = 1e5;
77
double Mu = E / (2 * (1 + nu)), Lam = E * nu / ((1 + nu) * (1 - 2 * nu));
88
double rho = 1000,
9-
k = 4e4, initial_stretch = 1, n_seg = 10, h = 0.01, side_len = 1, tol = 0.01, mu = 0.11;
9+
k = 4e4, initial_stretch = 1, n_seg = 4, h = 0.01, side_len = 1, tol = 0.01, mu = 0.11;
1010
// 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);
1111
InvFreeSimulator<double, 2> simulator(rho, side_len, initial_stretch, k, h, tol, mu, Mu, Lam, n_seg);
1212
simulator.run();

‎simulators/6_inv_free/src/simulator.cu‎

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -65,10 +65,10 @@ InvFreeSimulator<T, dim>::Impl::Impl(T rho, T side_len, T initial_stretch, T K,
6565
DBC.push_back((n_seg + 1) * (n_seg + 1));
6666
DBC_target.resize(DBC.size() * dim);
6767
DBC_limit.push_back(0);
68-
DBC_limit.push_back(-0.6);
68+
DBC_limit.push_back(-0.7);
6969
DBC_v.push_back(0);
70-
DBC_v.push_back(-0.3);
71-
DBC_stiff = 10;
70+
DBC_v.push_back(-0.5);
71+
DBC_stiff = 1000;
7272
x.push_back(0);
7373
x.push_back(side_len * 0.6);
7474
DBC_satisfied.resize(x.size() / dim);

‎simulators/7_self_contact/src/BarrierEnergy.cu‎

Lines changed: 61 additions & 30 deletions
Original file line numberDiff line numberDiff line change
@@ -14,6 +14,8 @@ struct BarrierEnergy<T, dim>::Impl
1414
DeviceBuffer<T> device_x;
1515
DeviceBuffer<T> device_contact_area, device_n, device_n_ceil, device_o;
1616
DeviceBuffer<int> device_bp, device_be;
17+
std::vector<int> bp;
18+
std::vector<int> be;
1719
int N;
1820
DeviceBuffer<T> device_grad;
1921
DeviceTripletMatrix<T, 1> device_hess;
@@ -48,6 +50,8 @@ BarrierEnergy<T, dim>::BarrierEnergy(const std::vector<T> &x, const std::vector<
4850
pimpl_->device_o.copy_from(o);
4951
pimpl_->device_bp.copy_from(bp);
5052
pimpl_->device_be.copy_from(be);
53+
pimpl_->bp = bp;
54+
pimpl_->be = be;
5155
pimpl_->device_hess.resize_triplets(pimpl_->N * dim * dim + (pimpl_->N - 1) * dim * dim * 4 + bp.size() * be.size() / 2 * 36);
5256
pimpl_->device_hess.reshape(x.size(), x.size());
5357
pimpl_->device_grad.resize(pimpl_->N * dim);
@@ -107,13 +111,14 @@ T BarrierEnergy<T, dim>::val()
107111
T d_sqr=PointEdgeDistanceVal(p,e0,e1);
108112
if(d_sqr<dhatsqr){
109113
T s = d_sqr / dhatsqr;
110-
device_val3(i)= kappa * device_contact_area(xI) * dhat/8*(s-1)*log(s);
114+
device_val3(i)= 0.5*kappa * device_contact_area(xI) * dhat/8*(s-1)*log(s);
111115
}
112116
} })
113117
.wait();
114-
return devicesum(device_val1) +
115-
devicesum(device_val2) +
116-
devicesum(device_val3);
118+
T val1 = devicesum(device_val1) +
119+
devicesum(device_val2);
120+
T val2 = devicesum(device_val3);
121+
return val1 + val2;
117122
} // Calculate the energy
118123

119124
template <typename T, int dim>
@@ -160,7 +165,7 @@ const DeviceBuffer<T> &BarrierEnergy<T, dim>::grad()
160165
for (int j = 0; j < dim; j++)
161166
{
162167
T grad =device_contact_area(i) * dhat * (kappa / 2 * (log(s) / dhat + (s - 1) / d)) * device_n_ceil(j);
163-
atomicAdd(&device_grad(i * dim + j), grad);
168+
device_grad(i * dim + j) += grad;
164169
atomicAdd(&device_grad((N-1) * dim + j), -grad);
165170
}
166171
} })
@@ -188,6 +193,7 @@ const DeviceBuffer<T> &BarrierEnergy<T, dim>::grad()
188193
}
189194
} })
190195
.wait();
196+
191197
return device_grad;
192198
}
193199

@@ -340,7 +346,6 @@ T BarrierEnergy<T, dim>::init_step_size(const DeviceBuffer<T> &p)
340346
alpha += device_n(j) * (device_x(i * dim + j) - device_o(j));
341347
}
342348
device_alpha(i) = min(device_alpha(i), 0.9 * alpha / -p_n);
343-
//printf("alpha: %f\n", device_alpha(i));
344349
} })
345350
.wait();
346351

@@ -351,7 +356,7 @@ T BarrierEnergy<T, dim>::init_step_size(const DeviceBuffer<T> &p)
351356
T p_n = 0;
352357
for (int j = 0; j < dim; j++)
353358
{
354-
p_n += p(i * dim + j) * device_n_ceil(j);
359+
p_n += (p(i * dim + j)-p((N-1)*dim+j)) * device_n_ceil(j);
355360
}
356361
if (p_n < 0)
357362
{
@@ -361,33 +366,59 @@ T BarrierEnergy<T, dim>::init_step_size(const DeviceBuffer<T> &p)
361366
alpha += device_n_ceil(j) * (device_x(i * dim + j) - device_x((N-1) * dim + j));
362367
}
363368
device_alpha(i) = min(device_alpha(i), 0.9 * alpha / -p_n);
364-
//printf("alpha: %f\n", device_alpha(i));
365369
} })
366370
.wait();
367371
T current_alpha = min_vector(device_alpha);
368-
ParallelFor(256)
369-
.apply(Npe, [current_alpha, device_x = device_x.cviewer(), P = p.cviewer(), device_alpha1 = device_alpha1.viewer(), device_bp = pimpl_->device_bp.cviewer(), device_be = pimpl_->device_be.cviewer(), Nbp, Nbe] __device__(int i) mutable
370-
{
371-
int xI = device_bp(i / Nbe);
372-
int eI0 = device_be(2*(i % Nbe)), eI1 = device_be(2*(i % Nbe) + 1);
373-
if (xI != eI0 && xI != eI1){
374-
Eigen::Matrix<T, 2, 1> p, e0, e1, dp, de0, de1;
375-
p<<device_x(xI*dim),device_x(xI*dim+1);
376-
e0<<device_x(eI0*dim),device_x(eI0*dim+1);
377-
e1<<device_x(eI1*dim),device_x(eI1*dim+1);
378-
dp<<P(xI*dim),P(xI*dim+1);
379-
de0<<P(eI0*dim),P(eI0*dim+1);
380-
de1<<P(eI1*dim),P(eI1*dim+1);
381-
if (bbox_overlap(p, e0, e1, dp, de0, de1, current_alpha))
382-
{
383-
T toc = narrow_phase_CCD(p, e0, e1, dp, de0, de1, current_alpha);
384-
printf("toc: %f\n", toc);
385-
device_alpha1(i) = min(device_alpha1(i), toc);
386-
}
387-
} })
388-
.wait();
372+
std::vector<T> x(device_x.size());
373+
device_x.copy_to(x);
374+
std::vector<T> cpu_p(p.size());
375+
p.copy_to(cpu_p);
376+
for (int i = 0; i < Nbp; i++)
377+
{
378+
for (int j = 0; j < Nbe; j++)
379+
{
380+
int xI = pimpl_->bp[i];
381+
int eI0 = pimpl_->be[2 * j], eI1 = pimpl_->be[2 * j + 1];
382+
if (xI != eI0 && xI != eI1)
383+
{
384+
Eigen::Matrix<T, 2, 1> p, e0, e1, dp, de0, de1;
385+
p << x[xI * dim], x[xI * dim + 1];
386+
e0 << x[eI0 * dim], x[eI0 * dim + 1];
387+
e1 << x[eI1 * dim], x[eI1 * dim + 1];
388+
dp << cpu_p[xI * dim], cpu_p[xI * dim + 1];
389+
de0 << cpu_p[eI0 * dim], cpu_p[eI0 * dim + 1];
390+
de1 << cpu_p[eI1 * dim], cpu_p[eI1 * dim + 1];
391+
if (bbox_overlap(p, e0, e1, dp, de0, de1, current_alpha))
392+
{
393+
T toc = narrow_phase_CCD(p, e0, e1, dp, de0, de1, current_alpha);
394+
if (toc < current_alpha)
395+
current_alpha = toc;
396+
}
397+
}
398+
}
399+
}
400+
// ParallelFor(256)
401+
// .apply(Npe, [current_alpha, device_x = device_x.cviewer(), P = p.cviewer(), device_alpha1 = device_alpha1.viewer(), device_bp = pimpl_->device_bp.cviewer(), device_be = pimpl_->device_be.cviewer(), Nbp, Nbe] __device__(int i) mutable
402+
// {
403+
// int xI = device_bp(i / Nbe);
404+
// int eI0 = device_be(2*(i % Nbe)), eI1 = device_be(2*(i % Nbe) + 1);
405+
// if (xI != eI0 && xI != eI1){
406+
// Eigen::Matrix<T, 2, 1> p, e0, e1, dp, de0, de1;
407+
// p<<device_x(xI*dim),device_x(xI*dim+1);
408+
// e0<<device_x(eI0*dim),device_x(eI0*dim+1);
409+
// e1<<device_x(eI1*dim),device_x(eI1*dim+1);
410+
// dp<<P(xI*dim),P(xI*dim+1);
411+
// de0<<P(eI0*dim),P(eI0*dim+1);
412+
// de1<<P(eI1*dim),P(eI1*dim+1);
413+
// if (bbox_overlap(p, e0, e1, dp, de0, de1, current_alpha))
414+
// {
415+
// T toc = narrow_phase_CCD(p, e0, e1, dp, de0, de1, current_alpha);
416+
// if (toc < device_alpha1(i))
417+
// device_alpha1(i) = toc;
418+
// }
419+
// } })
420+
// .wait();
389421
T a = min(min_vector(device_alpha1), current_alpha);
390-
printf("alpha: %f\n", a);
391422
return a;
392423
}
393424
template class BarrierEnergy<float, 2>;

‎simulators/7_self_contact/src/distance/CCD.cu‎

Lines changed: 0 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -34,10 +34,6 @@ __device__ __host__ T narrow_phase_CCD(Eigen::Matrix<T, 2, 1> p, Eigen::Matrix<T
3434
{
3535
return toc_upperbound;
3636
}
37-
printf("maxDispMag: %f\n", maxDispMag);
38-
printf("dp: %f, %f\n", dp[0], dp[1]);
39-
printf("de0: %f, %f\n", de0[0], de0[1]);
40-
printf("de1: %f, %f\n", de1[0], de1[1]);
4137
T eta = 0.1; // calculate the toc that first brings the distance to 0.1x the current distance
4238
T dist2_cur = PointEdgeDistanceVal(p, e0, e1);
4339
T dist_cur = std::sqrt(dist2_cur);
@@ -54,7 +50,6 @@ __device__ __host__ T narrow_phase_CCD(Eigen::Matrix<T, 2, 1> p, Eigen::Matrix<T
5450
e1 += tocLowerBound * de1;
5551
dist2_cur = PointEdgeDistanceVal(p, e0, e1);
5652
dist_cur = std::sqrt(dist2_cur);
57-
printf("dist_cur: %f\n", dist_cur);
5853
if (toc != 0 && dist_cur < gap)
5954
{
6055
break;

‎simulators/7_self_contact/src/main.cpp‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -6,7 +6,7 @@ int main()
66
double nu = 0.4, E = 1e5;
77
double Mu = E / (2 * (1 + nu)), Lam = E * nu / ((1 + nu) * (1 - 2 * nu));
88
double rho = 1000,
9-
k = 4e4, initial_stretch = 1, n_seg = 1, h = 0.01, side_len = 0.45, tol = 0.01, mu = 0.11;
9+
k = 4e4, initial_stretch = 1, n_seg = 2, h = 0.01, side_len = 0.45, tol = 0.01, mu = 0.11;
1010
// 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);
1111
InvFreeSimulator<double, 2> simulator(rho, side_len, initial_stretch, k, h, tol, mu, Mu, Lam, n_seg);
1212
simulator.run();

‎simulators/7_self_contact/src/simulator.cu‎

Lines changed: 9 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -64,7 +64,7 @@ InvFreeSimulator<T, dim>::Impl::Impl(T rho, T side_len, T initial_stretch, T K,
6464
generate(side_len, n_seg, x, e);
6565
for (int i = 0; i < (n_seg + 1) * (n_seg + 1); i++)
6666
{
67-
x.push_back(x[i * dim] - side_len * 0.05);
67+
x.push_back(x[i * dim] + side_len * 0.1);
6868
x.push_back(x[i * dim + 1] - side_len * 1.1);
6969
}
7070
int esize = e.size();
@@ -77,12 +77,12 @@ InvFreeSimulator<T, dim>::Impl::Impl(T rho, T side_len, T initial_stretch, T K,
7777
DBC.push_back(x.size() / dim);
7878
DBC_target.resize(DBC.size() * dim);
7979
DBC_limit.push_back(0);
80-
DBC_limit.push_back(-1.5);
80+
DBC_limit.push_back(-0.7);
8181
DBC_v.push_back(0);
8282
DBC_v.push_back(-0.5);
83-
DBC_stiff = 10;
83+
DBC_stiff = 1000;
8484
x.push_back(0);
85-
x.push_back(side_len * 1.2);
85+
x.push_back(side_len * 0.6);
8686
DBC_satisfied.resize(x.size() / dim);
8787
is_DBC.resize(x.size() / dim, 0);
8888
is_DBC[x.size() / dim - 1] = 1;
@@ -94,7 +94,7 @@ InvFreeSimulator<T, dim>::Impl::Impl(T rho, T side_len, T initial_stretch, T K,
9494
for (int i = 0; i < dim; i++)
9595
ground_n[i] /= n_norm;
9696
std::vector<T> ground_o(dim);
97-
ground_o[0] = 0, ground_o[1] = -1.7;
97+
ground_o[0] = 0, ground_o[1] = -1.0;
9898
v.resize(x.size(), 0);
9999
k.resize(e.size() / 2, K);
100100
l2.resize(e.size() / 2);
@@ -167,6 +167,7 @@ void InvFreeSimulator<T, dim>::Impl::step_forward()
167167
// std::cout << "Initial residual " << residual << "\n";
168168
while (residual > tol || DBC_satisfied.back() != 1) // use last one for simplisity, should check all
169169
{
170+
std::cout << "Iteration " << iter << " residual " << residual << " E_last" << E_last << "\n";
170171
if (residual <= tol && DBC_satisfied.back() != 1)
171172
{
172173
update_DBC_stiff(DBC_stiff * 2);
@@ -185,7 +186,7 @@ void InvFreeSimulator<T, dim>::Impl::step_forward()
185186
}
186187
std::cout << "step size = " << alpha << "\n";
187188
E_last = IP_val();
188-
std::cout << "Iteration " << iter << " residual " << residual << " E_last" << E_last << "\n";
189+
189190
p = search_direction();
190191
residual = max_vector(p) / h;
191192
iter += 1;
@@ -329,6 +330,8 @@ DeviceBuffer<T> InvFreeSimulator<T, dim>::Impl::search_direction()
329330
dir.resize(x.size());
330331
DeviceBuffer<T> grad = IP_grad();
331332
DeviceTripletMatrix<T, 1> hess = IP_hess();
333+
std::vector<T> hess_val(hess.values().size());
334+
hess.values().copy_to(hess_val.data());
332335
// check whether each DBC is satisfied
333336
for (int i = 0; i < DBC_satisfied.size(); i++)
334337
{

0 commit comments

Comments
 (0)