Skip to content

Commit 254ca00

Browse files
RoushelfyRoushelfy
committed
5 complete
1 parent 842a8a5 commit 254ca00

4 files changed

Lines changed: 61 additions & 67 deletions

File tree

‎simulators/5_mov_dirichlet/src/BarrierEnergy.cu‎

Lines changed: 30 additions & 38 deletions
Original file line numberDiff line numberDiff line change
@@ -40,10 +40,11 @@ BarrierEnergy<T, dim>::BarrierEnergy(const std::vector<T> &x, const std::vector<
4040
pimpl_->device_contact_area.copy_from(contact_area);
4141
std::vector<T> n_ceil(dim);
4242
n_ceil[1] = -1;
43+
n_ceil[0] = 0;
4344
pimpl_->device_n_ceil.copy_from(n_ceil);
4445
pimpl_->device_n.copy_from(n);
4546
pimpl_->device_o.copy_from(o);
46-
pimpl_->device_hess.resize_triplets((pimpl_->N * 2 - 1) * dim * dim);
47+
pimpl_->device_hess.resize_triplets(pimpl_->N * dim * dim + (pimpl_->N - 1) * dim * dim * 4);
4748
pimpl_->device_hess.reshape(x.size(), x.size());
4849
pimpl_->device_grad.resize(pimpl_->N * dim);
4950
}
@@ -129,7 +130,7 @@ const DeviceBuffer<T> &BarrierEnergy<T, dim>::grad()
129130
{
130131
T grad =device_contact_area(i) * dhat * (kappa / 2 * (log(s) / dhat + (s - 1) / d)) * device_n_ceil(j);
131132
device_grad(i * dim + j) += grad;
132-
device_grad((N-1) * dim + j) -= grad;
133+
atomicAdd(&device_grad((N-1) * dim + j), -grad);
133134
}
134135
} })
135136
.wait();
@@ -156,26 +157,18 @@ const DeviceTripletMatrix<T, 1> &BarrierEnergy<T, dim>::hess()
156157
{
157158
d += device_n(j) * (device_x(i * dim + j) - device_o(j));
158159
}
159-
if (d < dhat)
160-
for (int j = 0; j < dim; j++)
161-
{
162-
for (int k = 0; k < dim; k++)
163-
{
164-
int idx = i * dim * dim + j * dim + k;
165-
device_hess_row_idx(idx) = i * dim + j;
166-
device_hess_col_idx(idx) = i * dim + k;
167-
device_hess_val(idx) = device_contact_area(i) * dhat * kappa / (2 * d * d * dhat) * (d + dhat) * device_n(j) * device_n(k);
168-
}
169-
}
170-
else
160+
171161
for (int j = 0; j < dim; j++)
172162
{
173163
for (int k = 0; k < dim; k++)
174164
{
175165
int idx = i * dim * dim + j * dim + k;
176166
device_hess_row_idx(idx) = i * dim + j;
177167
device_hess_col_idx(idx) = i * dim + k;
178-
device_hess_val(idx) = 0;
168+
if (d < dhat)
169+
device_hess_val(idx) = device_contact_area(i) * dhat * kappa / (2 * d * d * dhat) * (d + dhat) * device_n(j) * device_n(k);
170+
else
171+
device_hess_val(idx) = 0;
179172
}
180173
} })
181174
.wait();
@@ -184,30 +177,27 @@ const DeviceTripletMatrix<T, 1> &BarrierEnergy<T, dim>::hess()
184177
T d = 0;
185178
for (int j = 0; j < dim; j++)
186179
{
187-
d += device_n_ceil(j) * (device_x(i * dim + j) - device_x((N-1) * dim + j));
180+
d += device_n_ceil(j) * (device_x(i * dim + j) - device_x((N - 1) * dim + j));
188181
}
189-
if (d < dhat)
190-
for (int j = 0; j < dim; j++)
191-
{
192-
for (int k = 0; k < dim; k++)
193-
{
194-
int idx =N*dim*dim+ i * dim * dim + j * dim + k;
195-
device_hess_row_idx(idx) = (N-1) * dim + j;
196-
device_hess_col_idx(idx) = (N-1) * dim + k;
197-
device_hess_val(idx) = device_contact_area(i) * dhat * kappa / (2 * d * d * dhat) * (d + dhat) * device_n_ceil(j) * device_n_ceil(k);
198-
}
199-
}
200-
else
201-
for (int j = 0; j < dim; j++)
202-
{
203-
for (int k = 0; k < dim; k++)
204-
{
205-
int idx = N*dim*dim+i * dim * dim + j * dim + k;
206-
device_hess_row_idx(idx) = (N-1) * dim + j;
207-
device_hess_col_idx(idx) = (N-1) * dim + k;
208-
device_hess_val(idx) = 0;
209-
}
210-
} })
182+
int index[2] = {i, N - 1};
183+
for (int nI = 0; nI < 2; nI++)
184+
for (int nJ = 0; nJ < 2; nJ++)
185+
for (int j = 0; j < dim; j++)
186+
{
187+
for (int k = 0; k < dim; k++)
188+
{
189+
int idx = N * dim * dim + i * dim * dim * 4 + nI * dim * dim * 2 + nJ * dim * dim + j * dim + k;
190+
device_hess_row_idx(idx) = index[nI] * dim + j;
191+
device_hess_col_idx(idx) = index[nJ] * dim + k;
192+
if (d < dhat)
193+
if (nI == nJ)
194+
device_hess_val(idx) = device_contact_area(i) * dhat * kappa / (2 * d * d * dhat) * (d + dhat) * device_n_ceil(j) * device_n_ceil(k);
195+
else
196+
device_hess_val(idx) = -device_contact_area(i) * dhat * kappa / (2 * d * d * dhat) * (d + dhat) * device_n_ceil(j) * device_n_ceil(k);
197+
else
198+
device_hess_val(idx) = 0;
199+
}
200+
} })
211201
.wait();
212202
return device_hess;
213203

@@ -240,6 +230,7 @@ T BarrierEnergy<T, dim>::init_step_size(const DeviceBuffer<T> &p)
240230
alpha += device_n(j) * (device_x(i * dim + j) - device_o(j));
241231
}
242232
device_alpha(i) = min(device_alpha(i), 0.9 * alpha / -p_n);
233+
//printf("alpha: %f\n", device_alpha(i));
243234
} })
244235
.wait();
245236

@@ -260,6 +251,7 @@ T BarrierEnergy<T, dim>::init_step_size(const DeviceBuffer<T> &p)
260251
alpha += device_n_ceil(j) * (device_x(i * dim + j) - device_x((N-1) * dim + j));
261252
}
262253
device_alpha(i) = min(device_alpha(i), 0.9 * alpha / -p_n);
254+
//printf("alpha: %f\n", device_alpha(i));
263255
} })
264256
.wait();
265257
return min_vector(device_alpha);

‎simulators/5_mov_dirichlet/src/SpringEnergy.cu‎

Lines changed: 10 additions & 13 deletions
Original file line numberDiff line numberDiff line change
@@ -46,7 +46,7 @@ SpringEnergy<T, dim>::SpringEnergy(const std::vector<T> &x, const std::vector<T>
4646
pimpl_->device_DBC_target.resize(DBC.size() * dim);
4747
pimpl_->k = k;
4848
pimpl_->device_grad.resize(pimpl_->N * dim);
49-
pimpl_->device_hess.resize_triplets(pimpl_->N * dim * dim);
49+
pimpl_->device_hess.resize_triplets(DBC.size() * dim);
5050
pimpl_->device_hess.reshape(x.size(), x.size());
5151
}
5252

@@ -81,17 +81,16 @@ T SpringEnergy<T, dim>::val()
8181

8282
ParallelFor(256).apply(N, [device_val = device_val.viewer(), device_x = device_x.cviewer(), device_m = device_m.cviewer(), device_DBC = device_DBC.cviewer(), device_DBC_target = device_DBC_target.cviewer(), k] __device__(int i) mutable
8383
{
84-
85-
int idx = device_DBC(i);
86-
T d = 0;
87-
for (int j = 0; j < dim; ++j)
88-
{
89-
d += (device_x(idx * dim + j) - device_DBC_target(i * dim + j)) * (device_x(idx * dim + j) - device_DBC_target(i * dim + j));
90-
}
91-
device_val(i) = 0.5 * k *d; })
84+
int idx = device_DBC(i);
85+
T d = 0;
86+
for (int j = 0; j < dim; ++j)
87+
{
88+
d += (device_x(idx * dim + j) - device_DBC_target(i * dim + j)) * (device_x(idx * dim + j) - device_DBC_target(i * dim + j));
89+
}
90+
device_val(i) = 0.5 * k * device_m(idx) * d; })
9291
.wait();
93-
94-
return devicesum(device_val);
92+
T ans = devicesum(device_val);
93+
return ans;
9594
}
9695

9796
template <typename T, int dim>
@@ -108,12 +107,10 @@ const DeviceBuffer<T> &SpringEnergy<T, dim>::grad()
108107

109108
ParallelFor(256).apply(N, [device_x = device_x.cviewer(), device_m = device_m.cviewer(), device_DBC = device_DBC.cviewer(), device_DBC_target = device_DBC_target.cviewer(), device_grad = device_grad.viewer(), k] __device__(int i) mutable
110109
{
111-
112110
int idx = device_DBC(i);
113111
for (int j = 0; j < dim; ++j)
114112
{
115113
device_grad(idx * dim + j) = (device_x(idx * dim + j) - device_DBC_target(i * dim + j))* k * device_m(idx);
116-
printf("grad[%d] = %f\n", idx * dim + j, device_grad(idx * dim + j));
117114
} })
118115
.wait();
119116

‎simulators/5_mov_dirichlet/src/main.cpp‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -2,7 +2,7 @@
22

33
int main()
44
{
5-
float rho = 1000, k = 2e4, initial_stretch = 1, n_seg = 4, h = 0.01, side_len = 1, tol = 0.01, mu = 0.11;
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);
77
MovDirichletSimulator<float, 2> simulator(rho, side_len, initial_stretch, k, h, tol, mu, n_seg);
88
simulator.run();

‎simulators/5_mov_dirichlet/src/simulator.cu‎

Lines changed: 20 additions & 15 deletions
Original file line numberDiff line numberDiff line change
@@ -17,7 +17,7 @@ struct MovDirichletSimulator<T, dim>::Impl
1717
T h, rho, side_len, initial_stretch, m, tol, mu, DBC_stiff;
1818
int resolution = 900, scale = 200, offset = resolution / 2, radius = 5;
1919
std::vector<T> x, x_tilde, v, k, l2, DBC_limit, DBC_v, DBC_target;
20-
std::vector<int> e, DBC, DBC_satisfied;
20+
std::vector<int> e, DBC, DBC_satisfied, is_DBC;
2121
DeviceBuffer<int> device_DBC;
2222
DeviceBuffer<T> device_contact_area;
2323
sf::RenderWindow window;
@@ -67,10 +67,13 @@ MovDirichletSimulator<T, dim>::Impl::Impl(T rho, T side_len, T initial_stretch,
6767
DBC_limit.push_back(0);
6868
DBC_limit.push_back(-0.6);
6969
DBC_v.push_back(0);
70-
DBC_v.push_back(-1);
70+
DBC_v.push_back(-0.5);
7171
DBC_stiff = 10;
7272
x.push_back(0);
7373
x.push_back(side_len * 0.6);
74+
DBC_satisfied.resize(x.size() / dim);
75+
is_DBC.resize(x.size() / dim, 0);
76+
is_DBC[(n_seg + 1) * (n_seg + 1)] = 1;
7477
std::vector<T> contact_area(x.size() / dim, side_len / n_seg);
7578
std::vector<T> ground_n(dim);
7679
ground_n[0] = 0, ground_n[1] = 1;
@@ -106,7 +109,7 @@ MovDirichletSimulator<T, dim>::Impl::Impl(T rho, T side_len, T initial_stretch,
106109
springenergy = SpringEnergy<T, dim>(x, std::vector<T>(N, m), DBC, DBC_stiff);
107110
DeviceBuffer<T> x_device(x);
108111
update_x(x_device);
109-
device_DBC = DeviceBuffer<int>(DBC);
112+
device_DBC = DeviceBuffer<int>(is_DBC);
110113
device_contact_area = DeviceBuffer<T>(contact_area);
111114
}
112115
template <typename T, int dim>
@@ -194,6 +197,7 @@ void MovDirichletSimulator<T, dim>::Impl::update_x(const DeviceBuffer<T> &new_x)
194197
massspringenergy.update_x(new_x);
195198
gravityenergy.update_x(new_x);
196199
barrierenergy.update_x(new_x);
200+
springenergy.update_x(new_x);
197201
new_x.copy_to(x);
198202
}
199203
template <typename T, int dim>
@@ -276,17 +280,13 @@ T MovDirichletSimulator<T, dim>::Impl::IP_val()
276280
template <typename T, int dim>
277281
DeviceBuffer<T> MovDirichletSimulator<T, dim>::Impl::IP_grad()
278282
{
279-
return add_vector<T>(add_vector<T>(add_vector<T>(add_vector<T>(add_vector<T>(inertialenergy.grad(),
280-
massspringenergy.grad(), 1.0, h * h),
281-
gravityenergy.grad(), 1.0, h * h),
282-
barrierenergy.grad(), 1.0, h * h),
283-
frictionenergy.grad(), 1.0, h * h),
284-
springenergy.grad(), 1.0, 1.0);
285-
// return add_vector<T>(add_vector<T>(add_vector<T>(add_vector<T>(inertialenergy.grad(),
286-
// massspringenergy.grad(), 1.0, h * h),
287-
// gravityenergy.grad(), 1.0, h * h),
288-
// barrierenergy.grad(), 1.0, h * h),
289-
// frictionenergy.grad(), 1.0, h * h);
283+
return add_vector<T>(
284+
add_vector<T>(add_vector<T>(add_vector<T>(add_vector<T>(inertialenergy.grad(),
285+
massspringenergy.grad(), 1.0, h * h),
286+
gravityenergy.grad(), 1.0, h * h),
287+
barrierenergy.grad(), 1.0, h * h),
288+
frictionenergy.grad(), 1.0, h * h),
289+
springenergy.grad(), 1.0, 1.0);
290290
}
291291

292292
template <typename T, int dim>
@@ -311,14 +311,19 @@ DeviceBuffer<T> MovDirichletSimulator<T, dim>::Impl::search_direction()
311311
DeviceBuffer<T> grad = IP_grad();
312312
DeviceTripletMatrix<T, 1> hess = IP_hess();
313313
// check whether each DBC is satisfied
314-
DBC_satisfied.resize(x.size() / dim, 0);
314+
for (int i = 0; i < DBC_satisfied.size(); i++)
315+
{
316+
DBC_satisfied[i] = 0;
317+
}
318+
315319
for (int i = 0; i < DBC.size(); i++)
316320
{
317321
T diff = 0;
318322
for (int d = 0; d < dim; d++)
319323
{
320324
diff += (x[DBC[i] * dim + d] - DBC_target[i * dim + d]) * (x[DBC[i] * dim + d] - DBC_target[i * dim + d]);
321325
}
326+
diff = sqrt(diff);
322327
if (diff / h < tol)
323328
{
324329
DBC_satisfied[DBC[i]] = 1;

0 commit comments

Comments
 (0)