Skip to content

Commit 55223c8

Browse files
AhdhnAhdhn
committed
fix NH
1 parent 71af8bd commit 55223c8

File tree

7 files changed

+76
-81
lines changed

7 files changed

+76
-81
lines changed

‎apps/NeoHookean/barrier_energy.h‎

Lines changed: 24 additions & 31 deletions
Original file line numberDiff line numberDiff line change
@@ -9,16 +9,15 @@ template <typename ProblemT,
99
typename VAttrT,
1010
typename VAttrI,
1111
typename T = typename VAttrT::Type>
12-
void floor_barrier_energy(ProblemT& problem,
13-
VAttrT& contact_area,
14-
const VertexHandle dbc_vertex,
15-
const VAttrT& x,
16-
const T h, // time_step
17-
const VAttrI& is_dbc,
18-
const vec3<T>& ground_n,
19-
const vec3<T>& ground_o,
20-
const T dhat,
21-
const T kappa)
12+
void floor_barrier_energy(ProblemT& problem,
13+
VAttrT& contact_area,
14+
const VAttrT& x,
15+
const T h, // time_step
16+
const VAttrI& is_dbc,
17+
const vec3<T>& ground_n,
18+
const vec3<T>& ground_o,
19+
const T dhat,
20+
const T kappa)
2221
{
2322

2423
const T h_sq = h * h;
@@ -32,9 +31,7 @@ void floor_barrier_energy(ProblemT& problem,
3231
[=] __device__(const auto& vh, auto& obj) mutable {
3332
using ActiveT = ACTIVE_TYPE(vh);
3433

35-
const Eigen::Vector3<T> x_dbc = x.template to_eigen<3>(dbc_vertex);
36-
37-
const Eigen::Vector3<ActiveT> xi = iter_val<ActiveT, 3>(vh, obj);
34+
const Eigen::Vector3<ActiveT> xi = iter_val<ActiveT, 3>(vh, x);
3835

3936
ActiveT E;
4037

@@ -45,44 +42,42 @@ void floor_barrier_energy(ProblemT& problem,
4542
if (d < dhat) {
4643
ActiveT s = d / dhat;
4744

45+
if (s <= T(0)) {
46+
using PassiveT = PassiveType<ActiveT>;
47+
return ActiveT(std::numeric_limits<PassiveT>::max());
48+
}
49+
4850
E = h_sq * contact_area(vh) * dhat * T(0.5) * kappa *
4951
(s - 1) * log(s);
5052
}
51-
52-
53-
// ceiling
54-
// d = (xi - x_dbc).dot(normal);
55-
//
56-
// if (d < dhat) {
57-
// ActiveT s = d / dhat;
58-
//
59-
// E += h_sq * contact_area(vh) * dhat * T(0.5) * kappa *
60-
// (s - 1) * log(s);
61-
//}
6253
}
6354

6455
return E;
6556
});
6657
}
6758

68-
template <typename VAttrT, typename T = typename VAttrT::Type>
59+
template <typename VAttrT, typename VAttrB, typename T = typename VAttrT::Type>
6960
void add_contact(RXMeshStatic& rx,
7061
CandidatePairsVV& contact_pairs,
7162
const VertexHandle dbc_vertex,
63+
const VAttrB& is_dbc,
7264
const VAttrT& x,
7365
const T dhat)
7466
{
7567
contact_pairs.reset();
7668

7769
rx.for_each_vertex(DEVICE, [=] __device__(const VertexHandle& vh) mutable {
70+
if (vh == dbc_vertex || is_dbc(vh)) {
71+
return;
72+
}
7873
const Eigen::Vector3<T> x_dbc = x.template to_eigen<3>(dbc_vertex);
7974
const Eigen::Vector3<T> xi = x.template to_eigen<3>(vh);
8075
const Eigen::Vector3<T> normal(0.0, -1.0, 0.0);
8176

8277
T d = (xi - x_dbc).dot(normal);
8378

8479
if (d < dhat) {
85-
contact_pairs.insert(vh, dbc_vertex);
80+
contact_pairs.insert(vh, dbc_vertex);
8681
}
8782
});
8883
}
@@ -101,9 +96,6 @@ void ceiling_barrier_energy(ProblemT& problem,
10196
{
10297
const T h_sq = h * h;
10398

104-
const Eigen::Vector3<T> o(ground_o[0], ground_o[1], ground_o[2]);
105-
const Eigen::Vector3<T> n(ground_n[0], ground_n[1], ground_n[2]);
106-
10799
const Eigen::Vector3<T> normal(0.0, -1.0, 0.0);
108100

109101
problem.template add_term<true>([=] __device__(const auto& id,
@@ -115,11 +107,11 @@ void ceiling_barrier_energy(ProblemT& problem,
115107

116108
const VertexHandle c1 = iter[1];
117109

110+
//???
111+
const Eigen::Vector3<ActiveT> xi = iter_val<ActiveT, 3>(id, iter, x, 0);
118112

119113
const Eigen::Vector3<T> x_dbc = x.template to_eigen<3>(c1);
120114

121-
const Eigen::Vector3<ActiveT> xi = iter_val<ActiveT, 3>(c0, obj);
122-
123115

124116
// ceiling
125117
ActiveT d = (xi - x_dbc).dot(normal);
@@ -175,6 +167,7 @@ T barrier_step_size(RXMeshStatic& rx,
175167
}
176168

177169
// ceiling
170+
//TODO this should be generalized
178171
if (!is_dbc(vh)) {
179172
p_n = glm::dot(n, (pi - p_dbc));
180173
if (p_n < 0) {

‎apps/NeoHookean/draw.h‎

Lines changed: 4 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -6,7 +6,7 @@ using namespace rxmesh;
66

77
template <typename FunT, typename VAttrT>
88
void draw(RXMeshStatic& rx,
9-
VAttrT& x_tilde,
9+
VAttrT& x,
1010
VAttrT& velocity,
1111
FunT& step_forward,
1212
int& step)
@@ -19,11 +19,11 @@ void draw(RXMeshStatic& rx,
1919
auto ps_callback = [&]() mutable {
2020
auto step_and_update = [&]() {
2121
step_forward();
22-
x_tilde.move(DEVICE, HOST);
22+
x.move(DEVICE, HOST);
2323
velocity.move(DEVICE, HOST);
2424
auto vel = rx.get_polyscope_mesh()->addVertexVectorQuantity(
2525
"Velocity", velocity);
26-
rx.get_polyscope_mesh()->updateVertexPositions(x_tilde);
26+
rx.get_polyscope_mesh()->updateVertexPositions(x);
2727
};
2828
if (ImGui::Button("Step")) {
2929
step_and_update();
@@ -41,8 +41,7 @@ void draw(RXMeshStatic& rx,
4141

4242
ImGui::SameLine();
4343
if (ImGui::Button("Export")) {
44-
rx.export_obj("NeoHookean_" + std::to_string(step) + ".obj",
45-
x_tilde);
44+
rx.export_obj("NeoHookean_" + std::to_string(step) + ".obj", x);
4645
}
4746

4847
if (is_running) {

‎apps/NeoHookean/gravity_energy.h‎

Lines changed: 12 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -11,22 +11,24 @@ void gravity_energy(ProblemT& problem,
1111
const T h,
1212
const T mass)
1313
{
14+
// obj here is x
15+
1416
const Eigen::Vector3<T> g(0.0, -9.81, 0.0);
1517

1618
const T neg_mass_times_h_sq = -mass * h * h;
1719

18-
problem.template add_term<Op::V, true>([=] __device__(const auto& vh,
19-
auto& obj) mutable {
20-
using ActiveT = ACTIVE_TYPE(vh);
21-
ActiveT E;
20+
problem.template add_term<Op::V, true>(
21+
[=] __device__(const auto& vh, auto& obj) mutable {
22+
using ActiveT = ACTIVE_TYPE(vh);
23+
ActiveT E;
2224

23-
if (int(is_dbc(vh)) == 0) {
24-
Eigen::Vector3<ActiveT> xi = iter_val<ActiveT, 3>(vh, obj);
25+
if (int(is_dbc(vh)) == 0) {
26+
Eigen::Vector3<ActiveT> xi = iter_val<ActiveT, 3>(vh, x);
2527

26-
E = neg_mass_times_h_sq * xi.dot(g);
27-
}
28+
E = neg_mass_times_h_sq * xi.dot(g);
29+
}
2830

2931

30-
return E;
31-
});
32+
return E;
33+
});
3234
}

‎apps/NeoHookean/inertial_energy.h‎

Lines changed: 4 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -7,6 +7,7 @@ using namespace rxmesh;
77
template <typename ProblemT, typename VAttrT, typename VAttrI, typename T>
88
void inertial_energy(ProblemT& problem,
99
const VAttrT& x,
10+
const VAttrT& x_tilde,
1011
const VAttrI& is_dbc,
1112
const T mass)
1213
{
@@ -19,11 +20,11 @@ void inertial_energy(ProblemT& problem,
1920

2021
if (is_dbc(vh) == 0) {
2122

22-
Eigen::Vector3<ActiveT> x_tilda = iter_val<ActiveT, 3>(vh, obj);
23+
Eigen::Vector3<ActiveT> xx = iter_val<ActiveT, 3>(vh, x);
2324

24-
Eigen::Vector3<T> xx = x.to_eigen<3>(vh);
25+
Eigen::Vector3<T> xx_tilde = x_tilde.to_eigen<3>(vh);
2526

26-
Eigen::Vector3<ActiveT> l = xx - x_tilda;
27+
Eigen::Vector3<ActiveT> l = xx - xx_tilde;
2728

2829
E = half_mass * l.squaredNorm();
2930
}

‎apps/NeoHookean/neo_hookean.cu‎

Lines changed: 25 additions & 28 deletions
Original file line numberDiff line numberDiff line change
@@ -81,6 +81,7 @@ void neo_hookean(RXMeshStatic& rx, T dx)
8181
volume.reset(0, DEVICE);
8282

8383
auto is_dbc_satisfied = *rx.add_vertex_attribute<int>("DBCSatisfied", 1);
84+
is_dbc_satisfied.reset(0, DEVICE);
8485

8586
VertexReduceHandle<int> rh(is_dbc_satisfied);
8687

@@ -107,12 +108,12 @@ void neo_hookean(RXMeshStatic& rx, T dx)
107108

108109
NetwtonSolver newton_solver(problem, &solver);
109110

110-
auto x = *rx.get_input_vertex_coordinates();
111+
auto& x = *problem.objective;
112+
x.copy_from(*rx.get_input_vertex_coordinates(), DEVICE, DEVICE);
111113

112-
auto x_n = *rx.add_vertex_attribute_like("x_n", x);
114+
auto x_n = *rx.add_vertex_attribute_like("x_n", x);
115+
auto x_tilde = *rx.add_vertex_attribute_like("x_tilde", x);
113116

114-
auto& x_tilde = *problem.objective;
115-
x_tilde.copy_from(x, DEVICE, DEVICE);
116117

117118
// Initializations
118119
init_volume_inverse_b(rx, x, volume, inv_b);
@@ -130,7 +131,7 @@ void neo_hookean(RXMeshStatic& rx, T dx)
130131

131132

132133
typename ProblemT::DenseMatT alpha(
133-
rx, std::max(rx.get_num_vertices(), rx.get_num_faces()), DEVICE);
134+
rx, std::max(rx.get_num_vertices(), rx.get_num_faces()), 1, DEVICE);
134135

135136
#if USE_POLYSCOPE
136137
// add BC to polyscope
@@ -142,10 +143,10 @@ void neo_hookean(RXMeshStatic& rx, T dx)
142143
#endif
143144

144145
// add inertial energy term OK
145-
inertial_energy(problem, x, is_dbc, mass);
146+
inertial_energy(problem, x, x_tilde, is_dbc, mass);
146147

147148
// add spring energy term
148-
spring_energy(problem, dbc_target, is_dbc, mass, dbc_stiff);
149+
spring_energy(problem, x, dbc_target, is_dbc, mass, dbc_stiff);
149150

150151

151152
// add gravity energy OK
@@ -154,7 +155,6 @@ void neo_hookean(RXMeshStatic& rx, T dx)
154155
// add barrier energy
155156
floor_barrier_energy(problem,
156157
contact_area,
157-
v_dbc[0],
158158
x,
159159
time_step,
160160
is_dbc,
@@ -206,27 +206,28 @@ void neo_hookean(RXMeshStatic& rx, T dx)
206206
x_n.copy_from(x, DEVICE, DEVICE);
207207

208208
// compute mu * lambda for each node using x_n
209-
compute_mu_lambda(rx,
209+
/*compute_mu_lambda(rx,
210210
fricition_coef,
211211
dhat,
212212
kappa,
213213
ground_n,
214214
ground_o,
215215
x,
216216
contact_area,
217-
mu_lambda);
217+
mu_lambda);*/
218218

219219
// target position for each DBC in the current time step
220220
update_dbc(
221221
rx, is_dbc, x, v_dbc_vel, v_dbc_limit, time_step, dbc_target);
222222

223223
// evaluate energy
224-
add_contact(rx, problem.vv_pairs, v_dbc[0], x_tilde, dhat);
224+
add_contact(rx, problem.vv_pairs, v_dbc[0], is_dbc, x, dhat);
225+
problem.update_hessian();
225226
problem.eval_terms();
226227

227228
// DBC satisfied
228229
check_dbc_satisfied(
229-
rx, is_dbc_satisfied, x_tilde, is_dbc, dbc_target, time_step, tol);
230+
rx, is_dbc_satisfied, x, is_dbc, dbc_target, time_step, tol);
230231

231232
// how many DBC are satisfied
232233
int num_satisfied = rh.reduce(is_dbc_satisfied, cub::Sum(), 0);
@@ -239,7 +240,6 @@ void neo_hookean(RXMeshStatic& rx, T dx)
239240
// get newton direction
240241
newton_solver.compute_direction();
241242

242-
243243
// residual is abs_max(newton_dir)/ h
244244
T residual = newton_solver.dir.abs_max() / time_step;
245245

@@ -249,7 +249,6 @@ void neo_hookean(RXMeshStatic& rx, T dx)
249249
int iter = 0;
250250
while (residual > tol || num_satisfied != num_dbc_vertices) {
251251

252-
253252
if (residual <= tol && num_satisfied != num_dbc_vertices) {
254253
dbc_stiff.multiply(T(2));
255254
}
@@ -271,17 +270,16 @@ void neo_hookean(RXMeshStatic& rx, T dx)
271270
newton_solver.line_search(line_search_init_step, 0.5);
272271

273272
// evaluate energy
274-
add_contact(rx, problem.vv_pairs, v_dbc[0], x_tilde, dhat);
273+
add_contact(rx, problem.vv_pairs, v_dbc[0], is_dbc, x, dhat);
274+
problem.update_hessian();
275275
problem.eval_terms();
276276

277+
// T f = problem.get_current_loss();
278+
// RXMESH_INFO("Subsetp, Energy: {}", f);
279+
277280
// DBC satisfied
278-
check_dbc_satisfied(rx,
279-
is_dbc_satisfied,
280-
x_tilde,
281-
is_dbc,
282-
dbc_target,
283-
time_step,
284-
tol);
281+
check_dbc_satisfied(
282+
rx, is_dbc_satisfied, x, is_dbc, dbc_target, time_step, tol);
285283

286284
// how many DBC are satisfied
287285
num_satisfied = rh.reduce(is_dbc_satisfied, cub::Sum(), 0);
@@ -301,13 +299,12 @@ void neo_hookean(RXMeshStatic& rx, T dx)
301299
// update velocity
302300
rx.for_each_vertex(
303301
DEVICE,
304-
[x, x_tilde, velocity, inv_time_step = 1.0 / time_step] __device__(
302+
[x, x_n, velocity, inv_time_step = 1.0 / time_step] __device__(
305303
VertexHandle vh) mutable {
306304
for (int i = 0; i < 3; ++i) {
307-
velocity(vh, i) =
308-
inv_time_step * (x_tilde(vh, i) - x(vh, i));
305+
velocity(vh, i) = inv_time_step * (x(vh, i) - x_n(vh, i));
309306

310-
x(vh, i) = x_tilde(vh, i);
307+
// x(vh, i) = x_tilde(vh, i);
311308
}
312309
});
313310

@@ -316,7 +313,7 @@ void neo_hookean(RXMeshStatic& rx, T dx)
316313
};
317314

318315
#if USE_POLYSCOPE
319-
draw(rx, x_tilde, velocity, step_forward, steps);
316+
draw(rx, x, velocity, step_forward, steps);
320317
#else
321318
while (steps < 5) {
322319
step_forward();
@@ -344,7 +341,7 @@ void neo_hookean(RXMeshStatic& rx, T dx)
344341

345342
int main(int argc, char** argv)
346343
{
347-
Log::init(spdlog::level::info);
344+
rx_init(0, spdlog::level::info);
348345

349346

350347
std::vector<std::vector<T>> verts;

‎apps/NeoHookean/neo_hookean_energy.h‎

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -34,9 +34,9 @@ void neo_hookean_energy(ProblemT& problem,
3434
return ActiveT();
3535
}
3636

37-
Eigen::Vector3<ActiveT> x0 = iter_val<ActiveT, 3>(fh, iter, obj, 0);
38-
Eigen::Vector3<ActiveT> x1 = iter_val<ActiveT, 3>(fh, iter, obj, 1);
39-
Eigen::Vector3<ActiveT> x2 = iter_val<ActiveT, 3>(fh, iter, obj, 2);
37+
Eigen::Vector3<ActiveT> x0 = iter_val<ActiveT, 3>(fh, iter, x, 0);
38+
Eigen::Vector3<ActiveT> x1 = iter_val<ActiveT, 3>(fh, iter, x, 1);
39+
Eigen::Vector3<ActiveT> x2 = iter_val<ActiveT, 3>(fh, iter, x, 2);
4040

4141
Eigen::Vector3<ActiveT> e0 = x2 - x0;
4242
Eigen::Vector3<ActiveT> e1 = x1 - x0;

0 commit comments

Comments
 (0)