Loading [MathJax]/jax/input/TeX/config.js
PolyFEM
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
SurfaceTractionForms.cpp
Go to the documentation of this file.
2
6
9#include <polyfem/State.hpp>
10
11// #include <finitediff.hpp>
12
13using namespace polyfem::utils;
14
15namespace polyfem::solver
16{
17 namespace
18 {
19 template <typename T>
20 Eigen::Matrix<T, Eigen::Dynamic, 1, 0, 3, 1> compute_displaced_normal(const Eigen::MatrixXd &reference_normal, const Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic, 0, 3, 3> &grad_x, const Eigen::MatrixXd &grad_u_local)
21 {
22 Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic, 0, 3, 3> trafo = grad_x;
23 for (int i = 0; i < grad_x.rows(); ++i)
24 for (int j = 0; j < grad_x.cols(); ++j)
25 trafo(i, j) = trafo(i, j) + grad_u_local(i, j);
26
27 Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic, 0, 3, 3> trafo_inv = inverse(trafo);
28
29 Eigen::Matrix<T, Eigen::Dynamic, 1, 0, 3, 1> n(reference_normal.cols(), 1);
30 for (int d = 0; d < n.size(); ++d)
31 n(d) = T(0);
32
33 for (int i = 0; i < n.size(); ++i)
34 for (int j = 0; j < n.size(); ++j)
35 n(j) = n(j) + (reference_normal(i) * trafo_inv(i, j));
36 n = n / n.norm();
37
38 return n;
39 }
40
41 class QuadraticBarrier : public ipc::Barrier
42 {
43 public:
44 QuadraticBarrier(const double weight = 1) : weight_(weight) {}
45
46 double operator()(const double d, const double dhat) const override
47 {
48 if (d > dhat)
49 return 0;
50 else
51 return weight_ * (d - dhat) * (d - dhat);
52 }
53 double first_derivative(const double d, const double dhat) const override
54 {
55 if (d > dhat)
56 return 0;
57 else
58 return 2 * weight_ * (d - dhat);
59 }
60 double second_derivative(const double d, const double dhat) const override
61 {
62 if (d > dhat)
63 return 0;
64 else
65 return 2 * weight_;
66 }
67 double units(const double dhat) const override
68 {
69 return dhat * dhat;
70 }
71
72 private:
73 const double weight_;
74 };
75
76 void compute_collision_mesh_quantities(
77 const State &state,
78 const std::set<int> &boundary_ids,
79 const ipc::CollisionMesh &collision_mesh,
80 Eigen::MatrixXd &node_positions,
81 Eigen::MatrixXi &boundary_edges,
82 Eigen::MatrixXi &boundary_triangles,
83 Eigen::SparseMatrix<double> &displacement_map,
84 std::vector<bool> &is_on_surface,
85 Eigen::MatrixXi &can_collide_cache)
86 {
87 std::vector<Eigen::Triplet<double>> displacement_map_entries;
89 node_positions, boundary_edges, boundary_triangles, displacement_map_entries);
90
91 is_on_surface.resize(node_positions.rows(), false);
92
93 std::map<int, std::set<int>> boundary_ids_to_dof;
95 Eigen::MatrixXd points, uv, normals;
96 Eigen::VectorXd weights;
97 Eigen::VectorXi global_primitive_ids;
98 for (const auto &lb : state.total_local_boundary)
99 {
100 const int e = lb.element_id();
101 bool has_samples = utils::BoundarySampler::boundary_quadrature(lb, state.n_boundary_samples(), *state.mesh, false, uv, points, normals, weights, global_primitive_ids);
102
103 if (!has_samples)
104 continue;
105
106 const basis::ElementBases &bs = state.bases[e];
107 const basis::ElementBases &gbs = state.geom_bases()[e];
108
109 vals.compute(e, state.mesh->is_volume(), points, bs, gbs);
110
111 for (int i = 0; i < lb.size(); ++i)
112 {
113 const int primitive_global_id = lb.global_primitive_id(i);
114 const auto nodes = bs.local_nodes_for_primitive(primitive_global_id, *state.mesh);
115 const int boundary_id = state.mesh->get_boundary_id(primitive_global_id);
116
117 if (!std::count(boundary_ids.begin(), boundary_ids.end(), boundary_id))
118 continue;
119
120 for (long n = 0; n < nodes.size(); ++n)
121 {
123 is_on_surface[v.global[0].index] = true;
124 if (v.global[0].index >= node_positions.rows())
125 log_and_throw_adjoint_error("Error building collision mesh in ProxyContactForceForm!");
126 boundary_ids_to_dof[boundary_id].insert(v.global[0].index);
127 }
128 }
129 }
130
131 if (!displacement_map_entries.empty())
132 {
133 displacement_map.resize(node_positions.rows(), state.n_bases);
134 displacement_map.setFromTriplets(displacement_map_entries.begin(), displacement_map_entries.end());
135 }
136
137 // Fix boundary edges and boundary triangles to exclude vertices not on triangles
138 Eigen::MatrixXi boundary_edges_alt(0, 2), boundary_triangles_alt(0, 3);
139 {
140 for (int i = 0; i < boundary_edges.rows(); ++i)
141 {
142 bool on_surface = true;
143 for (int j = 0; j < boundary_edges.cols(); ++j)
144 on_surface &= is_on_surface[boundary_edges(i, j)];
145 if (on_surface)
146 {
147 boundary_edges_alt.conservativeResize(boundary_edges_alt.rows() + 1, 2);
148 boundary_edges_alt.row(boundary_edges_alt.rows() - 1) = boundary_edges.row(i);
149 }
150 }
151
152 if (state.mesh->is_volume())
153 {
154 for (int i = 0; i < boundary_triangles.rows(); ++i)
155 {
156 bool on_surface = true;
157 for (int j = 0; j < boundary_triangles.cols(); ++j)
158 on_surface &= is_on_surface[boundary_triangles(i, j)];
159 if (on_surface)
160 {
161 boundary_triangles_alt.conservativeResize(boundary_triangles_alt.rows() + 1, 3);
162 boundary_triangles_alt.row(boundary_triangles_alt.rows() - 1) = boundary_triangles.row(i);
163 }
164 }
165 }
166 else
167 boundary_triangles_alt.resize(0, 0);
168 }
169 boundary_edges = boundary_edges_alt;
170 boundary_triangles = boundary_triangles_alt;
171
172 can_collide_cache.resize(0, 0);
173 can_collide_cache.resize(collision_mesh.num_vertices(), collision_mesh.num_vertices());
174 for (int i = 0; i < can_collide_cache.rows(); ++i)
175 {
176 int dof_idx_i = collision_mesh.to_full_vertex_id(i);
177 if (!is_on_surface[dof_idx_i])
178 continue;
179 for (int j = 0; j < can_collide_cache.cols(); ++j)
180 {
181 int dof_idx_j = collision_mesh.to_full_vertex_id(j);
182 if (!is_on_surface[dof_idx_j])
183 continue;
184
185 bool collision_allowed = true;
186 for (const auto &id : boundary_ids)
187 if (boundary_ids_to_dof[id].count(dof_idx_i) && boundary_ids_to_dof[id].count(dof_idx_j))
188 collision_allowed = false;
189 can_collide_cache(i, j) = collision_allowed;
190 }
191 }
192 }
193
194 double compute_quantity(
195 const Eigen::MatrixXd &node_positions,
196 const Eigen::MatrixXi &boundary_edges,
197 const Eigen::MatrixXi &boundary_triangles,
198 const Eigen::SparseMatrix<double> &displacement_map,
199 const std::vector<bool> &is_on_surface,
200 const Eigen::MatrixXi &can_collide_cache,
201 const double dhat,
202 const double dmin,
203 std::function<ipc::NormalCollisions(const Eigen::MatrixXd &)> cs_func,
204 const Eigen::MatrixXd &u,
205 const ipc::BarrierPotential &barrier_potential)
206
207 {
208 ipc::CollisionMesh collision_mesh = ipc::CollisionMesh(is_on_surface,
209 std::vector<bool>(is_on_surface.size(), false),
210 node_positions,
211 boundary_edges,
212 boundary_triangles,
213 displacement_map);
214
215 collision_mesh.can_collide = [&can_collide_cache](size_t vi, size_t vj) {
216 // return true;
217 return (bool)can_collide_cache(vi, vj);
218 };
219
220 collision_mesh.init_area_jacobians();
221
222 Eigen::MatrixXd displaced_surface = collision_mesh.displace_vertices(utils::unflatten(u, collision_mesh.dim()));
223
224 ipc::NormalCollisions cs_ = cs_func(displaced_surface);
225 cs_.build(collision_mesh, displaced_surface, dhat, dmin, ipc::build_broad_phase(ipc::BroadPhaseMethod::HASH_GRID));
226
227 Eigen::MatrixXd forces = collision_mesh.to_full_dof(barrier_potential.gradient(cs_, collision_mesh, displaced_surface));
228
229 double sum = (forces.array() * forces.array()).sum();
230
231 return sum;
232 }
233
234 } // namespace
235
237
238 // IntegrableFunctional TractionNormForm::get_integral_functional() const
239 // {
240 // IntegrableFunctional j;
241
242 // const std::string formulation = state_.formulation();
243 // const int power = in_power_;
244
245 // if (formulation == "Laplacian")
246 // log_and_throw_error("TractionNormForm is not implemented for Laplacian!");
247
248 // j.set_j([formulation, power, &state = std::as_const(state_)](const Eigen::MatrixXd &local_pts, const Eigen::MatrixXd &pts, const Eigen::MatrixXd &u, const Eigen::MatrixXd &grad_u, const Eigen::MatrixXd &lambda, const Eigen::MatrixXd &mu, const Eigen::MatrixXd &reference_normals, const assembler::ElementAssemblyValues &vals, const json &params, Eigen::MatrixXd &val) {
249 // val.setZero(grad_u.rows(), 1);
250 // int el_id = params["elem"];
251 // const double dt = state.problem->is_time_dependent() ? state.args["time"]["dt"].get<double>() : 0;
252 // const double t = params["t"];
253
254 // Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, 0, 3, 3> grad_u_local, grad_x;
255 // Eigen::MatrixXd grad_u_q, stress, traction_force, grad_unused;
256 // Eigen::MatrixXd reference_normal, displaced_normal;
257 // for (int q = 0; q < grad_u.rows(); q++)
258 // {
259 // grad_x = vals.jac_it[q].inverse();
260 // vector2matrix(grad_u.row(q), grad_u_q);
261 // grad_u_local = grad_u_q * grad_x;
262
263 // reference_normal = reference_normals.row(q);
264 // displaced_normal = compute_displaced_normal(reference_normal, grad_x, grad_u_local).transpose();
265 // // state.assembler->compute_stress_grad_multiply_vect(el_id, local_pts.row(q), pts.row(q), grad_u_q, Eigen::MatrixXd::Zero(1, grad_u_q.cols()), stress, grad_unused);
266 // state.assembler->compute_stress_grad_multiply_vect(OptAssemblerData(t, dt, el_id, local_pts.row(q), pts.row(q), grad_u_q), Eigen::MatrixXd::Zero(1, grad_u_q.cols()), stress, grad_unused);
267 // traction_force = displaced_normal * stress;
268 // val(q) = pow(traction_force.squaredNorm(), power / 2.);
269 // }
270 // });
271
272 // auto dj_dgradx = [formulation, power, &state = std::as_const(state_)](const Eigen::MatrixXd &local_pts, const Eigen::MatrixXd &pts, const Eigen::MatrixXd &u, const Eigen::MatrixXd &grad_u, const Eigen::MatrixXd &lambda, const Eigen::MatrixXd &mu, const Eigen::MatrixXd &reference_normals, const assembler::ElementAssemblyValues &vals, const json &params, Eigen::MatrixXd &val) {
273 // val.setZero(grad_u.rows(), grad_u.cols());
274 // int el_id = params["elem"];
275 // const double dt = state.problem->is_time_dependent() ? state.args["time"]["dt"].get<double>() : 0;
276 // const double t = params["t"];
277 // const int dim = sqrt(grad_u.cols());
278
279 // Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, 0, 3, 3> grad_u_local, grad_x;
280 // Eigen::MatrixXd grad_u_q, stress, traction_force, grad_unused;
281 // Eigen::MatrixXd reference_normal, displaced_normal;
282 // for (int q = 0; q < grad_u.rows(); q++)
283 // {
284 // grad_x = vals.jac_it[q].inverse();
285 // vector2matrix(grad_u.row(q), grad_u_q);
286 // grad_u_local = grad_u_q * grad_x;
287
288 // DiffScalarBase::setVariableCount(dim * dim);
289 // Eigen::Matrix<Diff, Eigen::Dynamic, Eigen::Dynamic, 0, 3, 3> grad_x_auto(dim, dim);
290 // for (int i = 0; i < dim; i++)
291 // for (int j = 0; j < dim; j++)
292 // grad_x_auto(i, j) = Diff(i + j * dim, grad_x(i, j));
293
294 // reference_normal = reference_normals.row(q);
295 // auto n = compute_displaced_normal(reference_normal, grad_x_auto, grad_u_local);
296 // displaced_normal.resize(1, dim);
297 // for (int i = 0; i < dim; ++i)
298 // displaced_normal(i) = n(i).getValue();
299 // // state.assembler->compute_stress_grad_multiply_vect(el_id, local_pts.row(q), pts.row(q), grad_u_q, Eigen::MatrixXd::Zero(1, grad_u_q.cols()), stress, grad_unused);
300 // state.assembler->compute_stress_grad_multiply_vect(OptAssemblerData(t, dt, el_id, local_pts.row(q), pts.row(q), grad_u_q), Eigen::MatrixXd::Zero(1, grad_u_q.cols()), stress, grad_unused);
301 // traction_force = displaced_normal * stress;
302
303 // const double coef = power * pow(traction_force.squaredNorm(), power / 2. - 1.);
304 // for (int k = 0; k < dim; ++k)
305 // for (int l = 0; l < dim; ++l)
306 // {
307 // double sum_j = 0;
308 // for (int j = 0; j < dim; ++j)
309 // {
310 // double grad_mult_stress = 0;
311 // for (int i = 0; i < dim; ++i)
312 // grad_mult_stress *= n(i).getGradient()(k + l * dim) * stress(i, j);
313
314 // sum_j += traction_force(j) * grad_mult_stress;
315 // }
316
317 // val(q, k * dim + l) = coef * sum_j;
318 // }
319 // }
320 // };
321 // j.set_dj_dgradx(dj_dgradx);
322 // j.set_dj_dgradu_local(dj_dgradx);
323
324 // auto dj_dgradu = [formulation, power, &state = std::as_const(state_)](const Eigen::MatrixXd &local_pts, const Eigen::MatrixXd &pts, const Eigen::MatrixXd &u, const Eigen::MatrixXd &grad_u, const Eigen::MatrixXd &lambda, const Eigen::MatrixXd &mu, const Eigen::MatrixXd &reference_normals, const assembler::ElementAssemblyValues &vals, const json &params, Eigen::MatrixXd &val) {
325 // val.setZero(grad_u.rows(), grad_u.cols());
326 // int el_id = params["elem"];
327 // const double dt = state.problem->is_time_dependent() ? state.args["time"]["dt"].get<double>() : 0;
328 // const double t = params["t"];
329 // const int dim = sqrt(grad_u.cols());
330
331 // Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, 0, 3, 3> grad_u_local, grad_x;
332 // Eigen::MatrixXd grad_u_q, stress, traction_force, vect_mult_dstress;
333 // Eigen::MatrixXd reference_normal, displaced_normal;
334 // for (int q = 0; q < grad_u.rows(); q++)
335 // {
336 // grad_x = vals.jac_it[q].inverse();
337 // vector2matrix(grad_u.row(q), grad_u_q);
338 // grad_u_local = grad_u_q * grad_x;
339
340 // reference_normal = reference_normals.row(q);
341 // displaced_normal = compute_displaced_normal(reference_normal, grad_x, grad_u_local).transpose();
342 // // state.assembler->compute_stress_grad_multiply_vect(el_id, local_pts.row(q), pts.row(q), grad_u_q, displaced_normal, stress, vect_mult_dstress);
343 // state.assembler->compute_stress_grad_multiply_vect(OptAssemblerData(t, dt, el_id, local_pts.row(q), pts.row(q), grad_u_q), displaced_normal, stress, vect_mult_dstress);
344 // traction_force = displaced_normal * stress;
345
346 // const double coef = power * pow(traction_force.squaredNorm(), power / 2. - 1.);
347 // for (int k = 0; k < dim; ++k)
348 // for (int l = 0; l < dim; ++l)
349 // {
350 // double sum_j = 0;
351 // for (int j = 0; j < dim; ++j)
352 // sum_j += traction_force(j) * vect_mult_dstress(l * dim + k, j);
353
354 // val(q, k * dim + l) = coef * sum_j;
355 // }
356 // }
357 // };
358 // j.set_dj_dgradu(dj_dgradu);
359
360 // /*
361 // j.set_dj_dgradu([formulation, power, &state = std::as_const(state_)](const Eigen::MatrixXd &local_pts, const Eigen::MatrixXd &pts, const Eigen::MatrixXd &u, const Eigen::MatrixXd &grad_u, const Eigen::MatrixXd &lambda, const Eigen::MatrixXd &mu, const Eigen::MatrixXd &reference_normals, const assembler::ElementAssemblyValues &vals, const json &params, Eigen::MatrixXd &val) {
362 // val.setZero(grad_u.rows(), grad_u.cols());
363 // int el_id = params["elem"];
364 // const int dim = sqrt(grad_u.cols());
365
366 // Eigen::MatrixXd displaced_normals;
367 // compute_displaced_normals(reference_normals, vals, grad_u, displaced_normals);
368
369 // Eigen::MatrixXd grad_u_q, stress, traction_force, normal_dstress;
370 // for (int q = 0; q < grad_u.rows(); q++)
371 // {
372 // vector2matrix(grad_u.row(q), grad_u_q);
373 // state.assembler->compute_stress_grad_multiply_vect(el_id, local_pts.row(q), pts.row(q), grad_u_q, displaced_normals.row(q), stress, normal_dstress);
374 // traction_force = displaced_normals.row(q) * stress;
375
376 // const double coef = power * pow(stress.squaredNorm(), power / 2. - 1.);
377 // for (int i = 0; i < dim; i++)
378 // for (int l = 0; l < dim; l++)
379 // val(q, i * dim + l) = coef * (traction_force.array() * normal_dstress.row(i * dim + l).array()).sum();
380 // }
381 // });
382
383 // auto dj_du = [formulation, power, &state = std::as_const(state_)](const Eigen::MatrixXd &local_pts, const Eigen::MatrixXd &pts, const Eigen::MatrixXd &u, const Eigen::MatrixXd &grad_u, const Eigen::MatrixXd &lambda, const Eigen::MatrixXd &mu, const Eigen::MatrixXd &reference_normals, const assembler::ElementAssemblyValues &vals, const json &params, Eigen::MatrixXd &val) {
384 // val.setZero(u.rows(), u.cols());
385 // int el_id = params["elem"];
386 // const int dim = sqrt(grad_u.cols());
387
388 // Eigen::MatrixXd displaced_normals;
389 // std::vector<Eigen::MatrixXd> normal_jacobian;
390 // compute_displaced_normal_jacobian(reference_normals, vals, grad_u, displaced_normals, normal_jacobian);
391
392 // Eigen::MatrixXd grad_u_q, stress, normal_duT_stress, traction_force, grad_unused;
393 // for (int q = 0; q < u.rows(); q++)
394 // {
395 // vector2matrix(grad_u.row(q), grad_u_q);
396 // state.assembler->compute_stress_grad_multiply_mat(el_id, local_pts.row(q), pts.row(q), grad_u_q, Eigen::MatrixXd::Zero(grad_u_q.rows(), grad_u_q.cols()), stress, grad_unused);
397 // traction_force = displaced_normals.row(q) * stress;
398
399 // Eigen::MatrixXd normal_du = normal_jacobian[q]; // compute this
400 // normal_duT_stress = normal_du.transpose() * stress;
401
402 // const double coef = power * pow(stress.squaredNorm(), power / 2. - 1.);
403 // for (int i = 0; i < dim; i++)
404 // val(q, i) = coef * (traction_force.array() * normal_duT_stress.row(i).array()).sum();
405 // }
406 // };
407 // j.set_dj_du(dj_du);
408 // j.set_dj_dx(dj_du);
409 // */
410
411 // /*
412 // const int normal_dim = 0;
413
414 // auto normal = [normal_dim](const Eigen::MatrixXd &reference_normal, const Eigen::MatrixXd &grad_x, const Eigen::MatrixXd &grad_u_local) {
415 // Eigen::MatrixXd trafo = grad_x + grad_u_local;
416 // Eigen::MatrixXd n = reference_normal * trafo.inverse();
417 // n.normalize();
418
419 // return n(normal_dim);
420 // };
421
422 // j.set_j([formulation, power, normal, normal_dim](const Eigen::MatrixXd &local_pts, const Eigen::MatrixXd &pts, const Eigen::MatrixXd &u, const Eigen::MatrixXd &grad_u, const Eigen::MatrixXd &lambda, const Eigen::MatrixXd &mu, const Eigen::MatrixXd &reference_normals, const assembler::ElementAssemblyValues &vals, const json &params, Eigen::MatrixXd &val) {
423 // val.setZero(grad_u.rows(), 1);
424 // int el_id = params["elem"];
425
426 // Eigen::MatrixXd grad_x, grad_u_local;
427 // for (int q = 0; q < grad_u.rows(); q++)
428 // {
429 // grad_x = vals.jac_it[q].inverse();
430 // vector2matrix(grad_u.row(q), grad_u_local);
431 // grad_u_local = grad_u_local * grad_x;
432 // double v = normal(reference_normals.row(q), grad_x, grad_u_local);
433 // val(q) = pow(v, 2);
434 // }
435 // });
436
437 // auto dj_dgradx = [formulation, power, normal](const Eigen::MatrixXd &local_pts, const Eigen::MatrixXd &pts, const Eigen::MatrixXd &u, const Eigen::MatrixXd &grad_u, const Eigen::MatrixXd &lambda, const Eigen::MatrixXd &mu, const Eigen::MatrixXd &reference_normals, const assembler::ElementAssemblyValues &vals, const json &params, Eigen::MatrixXd &val) {
438 // val.setZero(grad_u.rows(), grad_u.cols());
439 // int el_id = params["elem"];
440 // const int dim = sqrt(grad_u.cols());
441
442 // Eigen::MatrixXd grad_x, grad_u_local;
443 // for (int q = 0; q < grad_u.rows(); q++)
444 // {
445 // grad_x = vals.jac_it[q].inverse();
446 // vector2matrix(grad_u.row(q), grad_u_local);
447 // grad_u_local = grad_u_local * grad_x;
448
449 // const double v = normal(reference_normals.row(q), grad_x, grad_u_local);
450
451 // double eps = 1e-7;
452 // for (int i = 0; i < dim; ++i)
453 // for (int j = 0; j < dim; ++j)
454 // {
455 // Eigen::MatrixXd grad_x_copy = grad_x;
456 // grad_x_copy(i, j) += eps;
457 // const double val_plus = pow(normal(reference_normals.row(q), grad_x_copy, grad_u_local), 2);
458 // grad_x_copy(i, j) -= 2 * eps;
459 // const double val_minus = pow(normal(reference_normals.row(q), grad_x_copy, grad_u_local), 2);
460
461 // const double fd = (val_plus - val_minus) / (2 * eps);
462 // val(q, i * dim + j) = fd;
463 // }
464 // }
465 // };
466
467 // j.set_dj_dgradx(dj_dgradx);
468 // j.set_dj_dgradu_local(dj_dgradx);
469 // */
470
471 // // auto j_func = [formulation, power](const Eigen::MatrixXd &local_pts, const Eigen::MatrixXd &pts, const Eigen::MatrixXd &u, const Eigen::MatrixXd &grad_u, const Eigen::MatrixXd &lambda, const Eigen::MatrixXd &mu, const Eigen::MatrixXd &reference_normals, const assembler::ElementAssemblyValues &vals, const json &params, Eigen::MatrixXd &val) {
472 // // val.setZero(grad_u.rows(), 1);
473 // // int el_id = params["elem"];
474
475 // // Eigen::MatrixXd grad_u_q, stress, traction_force, grad_unused;
476 // // for (int q = 0; q < grad_u.rows(); q++)
477 // // {
478 // // vector2matrix(grad_u.row(q), grad_u_q);
479 // // double value = (grad_u_q * vals.jac_it[q].inverse() + vals.jac_it[q].inverse()).squaredNorm();
480 // // val(q) = pow(value, power / 2.);
481 // // }
482 // // };
483 // // j.set_j(j_func);
484
485 // // auto dj_dgradx = [formulation, power](const Eigen::MatrixXd &local_pts, const Eigen::MatrixXd &pts, const Eigen::MatrixXd &u, const Eigen::MatrixXd &grad_u, const Eigen::MatrixXd &lambda, const Eigen::MatrixXd &mu, const Eigen::MatrixXd &reference_normals, const assembler::ElementAssemblyValues &vals, const json &params, Eigen::MatrixXd &val) {
486 // // val.setZero(grad_u.rows(), grad_u.cols());
487 // // int el_id = params["elem"];
488 // // const int dim = sqrt(grad_u.cols());
489
490 // // Eigen::MatrixXd grad_u_q;
491 // // for (int q = 0; q < grad_u.rows(); q++)
492 // // {
493 // // vector2matrix(grad_u.row(q), grad_u_q);
494
495 // // Eigen::MatrixXd grad = (grad_u_q * vals.jac_it[q].inverse() + vals.jac_it[q].inverse());
496 // // const double coef = power * pow(grad.squaredNorm(), power / 2. - 1.);
497 // // val.row(q) = coef * flatten(grad);
498 // // }
499 // // };
500 // // j.set_dj_dgradx(dj_dgradx);
501 // // j.set_dj_dgradu_local(dj_dgradx);
502
503 // return j;
504 // }
505
506 // void TrueContactForceForm::build_active_nodes()
507 // {
508
509 // std::set<int> active_nodes_set = {};
510 // dim_ = state_.mesh->dimension();
511 // for (const auto &lb : state_.total_local_boundary)
512 // {
513 // const int e = lb.element_id();
514 // const basis::ElementBases &bs = state_.bases[e];
515
516 // for (int i = 0; i < lb.size(); i++)
517 // {
518 // const int global_primitive_id = lb.global_primitive_id(i);
519 // const auto nodes = bs.local_nodes_for_primitive(global_primitive_id, *state_.mesh);
520 // if (ids_.size() != 0 && ids_.find(state_.mesh->get_boundary_id(global_primitive_id)) == ids_.end())
521 // continue;
522
523 // for (long n = 0; n < nodes.size(); ++n)
524 // {
525 // const auto &b = bs.bases[nodes(n)];
526 // const int index = b.global()[0].index;
527
528 // for (int d = 0; d < dim_; ++d)
529 // active_nodes_set.insert(index * dim_ + d);
530 // }
531 // }
532 // }
533
534 // active_nodes_.resize(active_nodes_set.size());
535 // active_nodes_mat_.resize(state_.collision_mesh.full_ndof(), active_nodes_set.size());
536 // std::vector<Eigen::Triplet<double>> active_nodes_i;
537 // int count = 0;
538 // for (const auto node : active_nodes_set)
539 // {
540 // active_nodes_i.emplace_back(node, count, 1.0);
541 // active_nodes_(count++) = node;
542 // }
543 // active_nodes_mat_.setFromTriplets(active_nodes_i.begin(), active_nodes_i.end());
544
545 // epsv_ = state_.args["contact"]["epsv"];
546 // dhat_ = state_.args["contact"]["dhat"];
547 // friction_coefficient_ = state_.args["contact"]["friction_coefficient"];
548 // depends_on_step_prev_ = (friction_coefficient_ > 0);
549 // }
550
551 // double TrueContactForceForm::value_unweighted_step(const int time_step, const Eigen::VectorXd &x) const
552 // {
553 // assert(state_.solve_data.time_integrator != nullptr);
554 // assert(state_.solve_data.contact_form != nullptr);
555
556 // double barrier_stiffness = state_.solve_data.contact_form->weight();
557
558 // const ipc::CollisionMesh &collision_mesh = state_.collision_mesh;
559 // Eigen::MatrixXd displaced_surface = collision_mesh.displace_vertices(utils::unflatten(state_.diff_cached.u(time_step), dim_));
560
561 // Eigen::MatrixXd forces = Eigen::MatrixXd::Zero(collision_mesh.ndof(), 1);
562 // forces -= barrier_stiffness
563 // * state_.solve_data.contact_form->barrier_potential().gradient(
564 // state_.diff_cached.collision_set(time_step), collision_mesh, displaced_surface);
565 // if (state_.solve_data.friction_form && time_step > 0)
566 // {
567 // Eigen::MatrixXd surface_solution_prev = collision_mesh.vertices(utils::unflatten(state_.diff_cached.u(time_step - 1), collision_mesh.dim()));
568 // Eigen::MatrixXd surface_solution = collision_mesh.vertices(utils::unflatten(state_.diff_cached.u(time_step), collision_mesh.dim()));
569 // const Eigen::MatrixXd surface_velocities = (surface_solution - surface_solution_prev) / state_.solve_data.time_integrator->dt();
570
571 // forces += state_.solve_data.friction_form->friction_potential().force(
572 // state_.diff_cached.friction_collision_set(time_step),
573 // collision_mesh,
574 // collision_mesh.rest_positions(),
575 // surface_solution_prev,
576 // surface_velocities,
577 // state_.solve_data.contact_form->barrier_potential(),
578 // state_.solve_data.contact_form->barrier_stiffness(),
579 // state_.solve_data.friction_form->epsv());
580 // }
581 // forces = collision_mesh.to_full_dof(forces)(active_nodes_, Eigen::all);
582
583 // double sum = (forces.array() * forces.array()).sum();
584
585 // return sum;
586 // }
587
588 // Eigen::VectorXd TrueContactForceForm::compute_adjoint_rhs_unweighted_step(const int time_step, const Eigen::VectorXd &x, const State &state) const
589 // {
590 // assert(state_.solve_data.time_integrator != nullptr);
591 // assert(state_.solve_data.contact_form != nullptr);
592
593 // double barrier_stiffness = state_.solve_data.contact_form->weight();
594
595 // const ipc::CollisionMesh &collision_mesh = state_.collision_mesh;
596 // Eigen::MatrixXd displaced_surface = collision_mesh.displace_vertices(utils::unflatten(state_.diff_cached.u(time_step), dim_));
597
598 // Eigen::MatrixXd forces = Eigen::MatrixXd::Zero(collision_mesh.ndof(), 1);
599 // forces -= barrier_stiffness
600 // * state_.solve_data.contact_form->barrier_potential().gradient(
601 // state_.diff_cached.collision_set(time_step), collision_mesh, displaced_surface);
602 // if (state_.solve_data.friction_form && time_step > 0)
603 // {
604 // Eigen::MatrixXd surface_solution_prev = collision_mesh.vertices(utils::unflatten(state_.diff_cached.u(time_step - 1), collision_mesh.dim()));
605 // Eigen::MatrixXd surface_solution = collision_mesh.vertices(utils::unflatten(state_.diff_cached.u(time_step), collision_mesh.dim()));
606 // const Eigen::MatrixXd surface_velocities = (surface_solution - surface_solution_prev) / state_.solve_data.time_integrator->dt();
607
608 // forces += state_.solve_data.friction_form->friction_potential().force(
609 // state_.diff_cached.friction_collision_set(time_step),
610 // collision_mesh,
611 // collision_mesh.rest_positions(),
612 // surface_solution_prev,
613 // surface_velocities,
614 // state_.solve_data.contact_form->barrier_potential(),
615 // state_.solve_data.contact_form->barrier_stiffness(),
616 // state_.solve_data.friction_form->epsv());
617 // }
618 // forces = collision_mesh.to_full_dof(forces)(active_nodes_, Eigen::all);
619
620 // StiffnessMatrix hessian(collision_mesh.ndof(), collision_mesh.ndof());
621 // hessian -= barrier_stiffness
622 // * state_.solve_data.contact_form->barrier_potential().hessian(
623 // state_.diff_cached.collision_set(time_step), collision_mesh, displaced_surface, false);
624 // if (state_.solve_data.friction_form && time_step > 0)
625 // {
626 // Eigen::MatrixXd surface_solution_prev = collision_mesh.vertices(utils::unflatten(state_.diff_cached.u(time_step - 1), collision_mesh.dim()));
627 // Eigen::MatrixXd surface_solution = collision_mesh.vertices(utils::unflatten(state_.diff_cached.u(time_step), collision_mesh.dim()));
628 // const Eigen::MatrixXd surface_velocities = (surface_solution - surface_solution_prev) / state_.solve_data.time_integrator->dt();
629
630 // const double dv_du = -1 / state_.solve_data.time_integrator->dt();
631 // hessian += state_.solve_data.friction_form->friction_potential().force_jacobian(
632 // state_.diff_cached.friction_collision_set(time_step),
633 // collision_mesh,
634 // collision_mesh.rest_positions(),
635 // surface_solution_prev,
636 // surface_velocities,
637 // state_.solve_data.contact_form->barrier_potential(),
638 // state_.solve_data.contact_form->barrier_stiffness(),
639 // ipc::FrictionPotential::DiffWRT::VELOCITIES);
640 // }
641 // hessian = collision_mesh.to_full_dof(hessian);
642
643 // Eigen::VectorXd gradu = 2 * hessian.transpose() * active_nodes_mat_ * forces;
644
645 // return gradu;
646 // }
647
648 // Eigen::VectorXd TrueContactForceForm::compute_adjoint_rhs_unweighted_step_prev(const int time_step, const Eigen::VectorXd &x, const State &state) const
649 // {
650 // assert(state_.solve_data.time_integrator != nullptr);
651 // assert(state_.solve_data.contact_form != nullptr);
652
653 // double barrier_stiffness = state_.solve_data.contact_form->weight();
654
655 // const ipc::CollisionMesh &collision_mesh = state_.collision_mesh;
656 // Eigen::MatrixXd displaced_surface = collision_mesh.displace_vertices(utils::unflatten(state_.diff_cached.u(time_step), dim_));
657
658 // StiffnessMatrix hessian_prev(collision_mesh.ndof(), collision_mesh.ndof());
659 // Eigen::MatrixXd forces = Eigen::MatrixXd::Zero(collision_mesh.ndof(), 1);
660
661 // if (state_.solve_data.friction_form && time_step > 0)
662 // {
663 // Eigen::MatrixXd surface_solution_prev = collision_mesh.vertices(utils::unflatten(state_.diff_cached.u(time_step - 1), collision_mesh.dim()));
664 // Eigen::MatrixXd surface_solution = collision_mesh.vertices(utils::unflatten(state_.diff_cached.u(time_step), collision_mesh.dim()));
665 // const Eigen::MatrixXd surface_velocities = (surface_solution - surface_solution_prev) / state_.solve_data.time_integrator->dt();
666
667 // forces -= barrier_stiffness
668 // * state_.solve_data.contact_form->barrier_potential().gradient(
669 // state_.diff_cached.collision_set(time_step), collision_mesh, displaced_surface);
670 // if (state_.solve_data.friction_form && time_step > 0)
671 // forces += state_.solve_data.friction_form->friction_potential().force(
672 // state_.diff_cached.friction_collision_set(time_step),
673 // collision_mesh,
674 // collision_mesh.rest_positions(),
675 // surface_solution_prev,
676 // surface_velocities,
677 // state_.solve_data.contact_form->barrier_potential(),
678 // state_.solve_data.contact_form->barrier_stiffness(),
679 // state_.solve_data.friction_form->epsv());
680 // forces = collision_mesh.to_full_dof(forces)(active_nodes_, Eigen::all);
681
682 // const double dv_du = -1 / state_.solve_data.time_integrator->dt();
683 // hessian_prev += state_.solve_data.friction_form->friction_potential().force_jacobian(
684 // state_.diff_cached.friction_collision_set(time_step),
685 // collision_mesh,
686 // collision_mesh.rest_positions(),
687 // surface_solution_prev,
688 // surface_velocities,
689 // state_.solve_data.contact_form->barrier_potential(),
690 // state_.solve_data.contact_form->barrier_stiffness(),
691 // ipc::FrictionPotential::DiffWRT::LAGGED_DISPLACEMENTS);
692 // hessian_prev += dv_du
693 // * state_.solve_data.friction_form->friction_potential().force_jacobian(
694 // state_.diff_cached.friction_collision_set(time_step),
695 // collision_mesh,
696 // collision_mesh.rest_positions(),
697 // surface_solution_prev,
698 // surface_velocities,
699 // state_.solve_data.contact_form->barrier_potential(),
700 // state_.solve_data.contact_form->barrier_stiffness(),
701 // ipc::FrictionPotential::DiffWRT::VELOCITIES);
702 // }
703 // hessian_prev = collision_mesh.to_full_dof(hessian_prev);
704
705 // Eigen::VectorXd gradu = 2 * hessian_prev.transpose() * active_nodes_mat_ * forces;
706
707 // return gradu;
708 // }
709
710 // void TrueContactForceForm::compute_partial_gradient_unweighted_step(const int time_step, const Eigen::VectorXd &x, Eigen::VectorXd &gradv) const
711 // {
712 // assert(state_.solve_data.time_integrator != nullptr);
713 // assert(state_.solve_data.contact_form != nullptr);
714
715 // double barrier_stiffness = state_.solve_data.contact_form->weight();
716
717 // const ipc::CollisionMesh &collision_mesh = state_.collision_mesh;
718 // Eigen::MatrixXd displaced_surface = collision_mesh.displace_vertices(utils::unflatten(state_.diff_cached.u(time_step), dim_));
719
720 // Eigen::MatrixXd forces = Eigen::MatrixXd::Zero(collision_mesh.ndof(), 1);
721 // forces -= barrier_stiffness
722 // * state_.solve_data.contact_form->barrier_potential().gradient(
723 // state_.diff_cached.collision_set(time_step), collision_mesh, displaced_surface);
724 // if (state_.solve_data.friction_form && time_step > 0)
725 // {
726 // Eigen::MatrixXd surface_solution_prev = collision_mesh.vertices(utils::unflatten(state_.diff_cached.u(time_step - 1), collision_mesh.dim()));
727 // Eigen::MatrixXd surface_solution = collision_mesh.vertices(utils::unflatten(state_.diff_cached.u(time_step), collision_mesh.dim()));
728 // const Eigen::MatrixXd surface_velocities = (surface_solution - surface_solution_prev) / state_.solve_data.time_integrator->dt();
729
730 // forces += state_.solve_data.friction_form->friction_potential().force(
731 // state_.diff_cached.friction_collision_set(time_step),
732 // collision_mesh,
733 // collision_mesh.rest_positions(),
734 // surface_solution_prev,
735 // surface_velocities,
736 // state_.solve_data.contact_form->barrier_potential(),
737 // state_.solve_data.contact_form->barrier_stiffness(),
738 // state_.solve_data.friction_form->epsv());
739 // }
740 // forces = collision_mesh.to_full_dof(forces)(active_nodes_, Eigen::all);
741
742 // StiffnessMatrix shape_derivative(collision_mesh.ndof(), collision_mesh.ndof());
743 // shape_derivative -= barrier_stiffness
744 // * state_.solve_data.contact_form->barrier_potential().shape_derivative(
745 // state_.diff_cached.collision_set(time_step), collision_mesh, displaced_surface);
746 // if (state_.solve_data.friction_form && time_step > 0)
747 // {
748 // Eigen::MatrixXd surface_solution_prev = collision_mesh.vertices(utils::unflatten(state_.diff_cached.u(time_step - 1), collision_mesh.dim()));
749 // Eigen::MatrixXd surface_solution = collision_mesh.vertices(utils::unflatten(state_.diff_cached.u(time_step), collision_mesh.dim()));
750 // const Eigen::MatrixXd surface_velocities = (surface_solution - surface_solution_prev) / state_.solve_data.time_integrator->dt();
751
752 // shape_derivative += state_.solve_data.friction_form->friction_potential().force_jacobian(
753 // state_.diff_cached.friction_collision_set(time_step),
754 // collision_mesh,
755 // collision_mesh.rest_positions(),
756 // surface_solution_prev,
757 // surface_velocities,
758 // state_.solve_data.contact_form->barrier_potential(),
759 // state_.solve_data.contact_form->barrier_stiffness(),
760 // ipc::FrictionPotential::DiffWRT::REST_POSITIONS);
761 // }
762 // shape_derivative = collision_mesh.to_full_dof(shape_derivative);
763
764 // Eigen::VectorXd grads = 2 * shape_derivative.transpose() * active_nodes_mat_ * forces;
765 // grads = state_.basis_nodes_to_gbasis_nodes * grads;
766 // grads = AdjointTools::map_node_to_primitive_order(state_, grads);
767
768 // gradv.setZero(x.size());
769
770 // for (const auto &param_map : variable_to_simulations_)
771 // {
772 // const auto &param_type = param_map->get_parameter_type();
773
774 // for (const auto &state : param_map->get_states())
775 // {
776 // if (state.get() != &state_)
777 // continue;
778
779 // if (param_type != ParameterType::Shape)
780 // log_and_throw_error("Only support contact force derivative wrt. shape!");
781
782 // if (grads.size() > 0)
783 // gradv += param_map->apply_parametrization_jacobian(grads, x);
784 // }
785 // }
786 // }
787
789 const VariableToSimulationGroup &variable_to_simulations,
790 const State &state,
791 const double dhat,
792 const bool quadratic_potential,
793 const json &args)
794 : StaticForm(variable_to_simulations),
795 state_(state),
796 dhat_(dhat),
797 dmin_(0),
798 barrier_potential_(dhat, true)
799 {
800 auto tmp_ids = args["surface_selection"].get<std::vector<int>>();
801 boundary_ids_ = std::set(tmp_ids.begin(), tmp_ids.end());
802
804
805 broad_phase_method_ = ipc::BroadPhaseMethod::HASH_GRID;
806
807 if (state.problem->is_time_dependent())
808 {
809 int time_steps = state.args["time"]["time_steps"].get<int>() + 1;
810 collision_set_indicator_.setZero(time_steps);
811 for (int i = 0; i < time_steps + 1; ++i)
812 {
813 collision_sets_.push_back(std::make_shared<ipc::NormalCollisions>());
814 // collision_sets_.back()->set_use_convergent_formulation(true);
815 collision_sets_.back()->set_use_area_weighting(true);
816 collision_sets_.back()->set_use_improved_max_approximator(true);
817 collision_sets_.back()->set_enable_shape_derivatives(true);
818 }
819 }
820 else
821 {
822 collision_set_indicator_.setZero(1);
823 collision_sets_.push_back(std::make_shared<ipc::NormalCollisions>());
824 collision_sets_.back()->set_use_area_weighting(true);
825 collision_sets_.back()->set_use_improved_max_approximator(true);
826 collision_sets_.back()->set_enable_shape_derivatives(true);
827 }
828
829 if (quadratic_potential)
830 barrier_potential_.set_barrier(std::make_shared<QuadraticBarrier>());
831 }
832
834 {
835 boundary_ids_to_dof_.clear();
836 can_collide_cache_.resize(0, 0);
837 node_positions_.resize(0, 0);
838
839 // Eigen::MatrixXd node_positions;
840 Eigen::MatrixXi boundary_edges, boundary_triangles;
841 std::vector<Eigen::Triplet<double>> displacement_map_entries;
843 node_positions_, boundary_edges, boundary_triangles, displacement_map_entries);
844
845 std::vector<bool> is_on_surface;
846 is_on_surface.resize(node_positions_.rows(), false);
847
849 Eigen::MatrixXd points, uv, normals;
850 Eigen::VectorXd weights;
851 Eigen::VectorXi global_primitive_ids;
852 for (const auto &lb : state_.total_local_boundary)
853 {
854 const int e = lb.element_id();
855 bool has_samples = utils::BoundarySampler::boundary_quadrature(lb, state_.n_boundary_samples(), *state_.mesh, false, uv, points, normals, weights, global_primitive_ids);
856
857 if (!has_samples)
858 continue;
859
860 const basis::ElementBases &bs = state_.bases[e];
861 const basis::ElementBases &gbs = state_.geom_bases()[e];
862
863 vals.compute(e, state_.mesh->is_volume(), points, bs, gbs);
864
865 for (int i = 0; i < lb.size(); ++i)
866 {
867 const int primitive_global_id = lb.global_primitive_id(i);
868 const auto nodes = bs.local_nodes_for_primitive(primitive_global_id, *state_.mesh);
869 const int boundary_id = state_.mesh->get_boundary_id(primitive_global_id);
870
871 if (!std::count(boundary_ids_.begin(), boundary_ids_.end(), boundary_id))
872 continue;
873
874 for (long n = 0; n < nodes.size(); ++n)
875 {
876 const assembler::AssemblyValues &v = vals.basis_values[nodes(n)];
877 is_on_surface[v.global[0].index] = true;
878 if (v.global[0].index >= node_positions_.rows())
879 log_and_throw_adjoint_error("Error building collision mesh in ProxyContactForceForm!");
880 boundary_ids_to_dof_[boundary_id].insert(v.global[0].index);
881 }
882 }
883 }
884
885 Eigen::SparseMatrix<double> displacement_map;
886 if (!displacement_map_entries.empty())
887 {
888 displacement_map.resize(node_positions_.rows(), state_.n_bases);
889 displacement_map.setFromTriplets(displacement_map_entries.begin(), displacement_map_entries.end());
890 }
891
892 // Fix boundary edges and boundary triangles to exclude vertices not on triangles
893 Eigen::MatrixXi boundary_edges_alt(0, 2), boundary_triangles_alt(0, 3);
894 {
895 for (int i = 0; i < boundary_edges.rows(); ++i)
896 {
897 bool on_surface = true;
898 for (int j = 0; j < boundary_edges.cols(); ++j)
899 on_surface &= is_on_surface[boundary_edges(i, j)];
900 if (on_surface)
901 {
902 boundary_edges_alt.conservativeResize(boundary_edges_alt.rows() + 1, 2);
903 boundary_edges_alt.row(boundary_edges_alt.rows() - 1) = boundary_edges.row(i);
904 }
905 }
906
907 if (state_.mesh->is_volume())
908 {
909 for (int i = 0; i < boundary_triangles.rows(); ++i)
910 {
911 bool on_surface = true;
912 for (int j = 0; j < boundary_triangles.cols(); ++j)
913 on_surface &= is_on_surface[boundary_triangles(i, j)];
914 if (on_surface)
915 {
916 boundary_triangles_alt.conservativeResize(boundary_triangles_alt.rows() + 1, 3);
917 boundary_triangles_alt.row(boundary_triangles_alt.rows() - 1) = boundary_triangles.row(i);
918 }
919 }
920 }
921 else
922 boundary_triangles_alt.resize(0, 0);
923 }
924
925 collision_mesh_ = ipc::CollisionMesh(is_on_surface,
926 std::vector<bool>(is_on_surface.size(), false),
928 boundary_edges_alt,
929 boundary_triangles_alt,
930 displacement_map);
931
932 can_collide_cache_.resize(collision_mesh_.num_vertices(), collision_mesh_.num_vertices());
933 for (int i = 0; i < can_collide_cache_.rows(); ++i)
934 {
935 int dof_idx_i = collision_mesh_.to_full_vertex_id(i);
936 if (!is_on_surface[dof_idx_i])
937 continue;
938 for (int j = 0; j < can_collide_cache_.cols(); ++j)
939 {
940 int dof_idx_j = collision_mesh_.to_full_vertex_id(j);
941 if (!is_on_surface[dof_idx_j])
942 continue;
943
944 bool collision_allowed = true;
945 for (const auto &id : boundary_ids_)
946 if (boundary_ids_to_dof_[id].count(dof_idx_i) && boundary_ids_to_dof_[id].count(dof_idx_j))
947 collision_allowed = false;
948 can_collide_cache_(i, j) = collision_allowed;
949 }
950 }
951
952 collision_mesh_.can_collide = [this](size_t vi, size_t vj) {
953 // return true;
954 return (bool)can_collide_cache_(vi, vj);
955 };
956
957 collision_mesh_.init_area_jacobians();
958 }
959
960 const ipc::NormalCollisions &ProxyContactForceForm::get_or_compute_collision_set(const int time_step, const Eigen::MatrixXd &displaced_surface) const
961 {
962 if (!collision_set_indicator_(time_step))
963 {
964 collision_sets_[time_step]->build(
965 collision_mesh_, displaced_surface, dhat_, dmin_, ipc::build_broad_phase(broad_phase_method_));
966 collision_set_indicator_(time_step) = 1;
967 }
968 return *collision_sets_[time_step];
969 }
970
971 double ProxyContactForceForm::value_unweighted_step(const int time_step, const Eigen::VectorXd &x) const
972 {
973 assert(state_.solve_data.time_integrator != nullptr);
974 assert(state_.solve_data.contact_form != nullptr);
975
976 const Eigen::MatrixXd displaced_surface = collision_mesh_.displace_vertices(utils::unflatten(state_.diff_cached.u(time_step), collision_mesh_.dim()));
977 auto collision_set = get_or_compute_collision_set(time_step, displaced_surface);
978
979 Eigen::MatrixXd forces = collision_mesh_.to_full_dof(barrier_potential_.gradient(collision_set, collision_mesh_, displaced_surface));
980
981 double sum = (forces.array() * forces.array()).sum();
982
983 return sum;
984 }
985
986 Eigen::VectorXd ProxyContactForceForm::compute_adjoint_rhs_step(const int time_step, const Eigen::VectorXd &x, const State &state) const
987 {
988 assert(state_.solve_data.time_integrator != nullptr);
989 assert(state_.solve_data.contact_form != nullptr);
990
991 const Eigen::MatrixXd displaced_surface = collision_mesh_.displace_vertices(utils::unflatten(state_.diff_cached.u(time_step), collision_mesh_.dim()));
992 auto collision_set = get_or_compute_collision_set(time_step, displaced_surface);
993
994 Eigen::MatrixXd forces = collision_mesh_.to_full_dof(barrier_potential_.gradient(collision_set, collision_mesh_, displaced_surface));
995 StiffnessMatrix hessian = collision_mesh_.to_full_dof(barrier_potential_.hessian(collision_set, collision_mesh_, displaced_surface, ipc::PSDProjectionMethod::NONE));
996
997 Eigen::VectorXd gradu = 2 * hessian.transpose() * forces;
998
999 // logger().trace("u norm {}", state_.diff_cached.u(time_step).norm());
1000
1001 // Eigen::VectorXd G;
1002 // fd::finite_gradient(
1003 // state_.diff_cached.u(time_step), [this, time_step](const Eigen::VectorXd &x_) {
1004 // const Eigen::MatrixXd displaced_surface = collision_mesh_.displace_vertices(utils::unflatten(x_, collision_mesh_.dim()));
1005
1006 // auto collision_set = get_or_compute_collision_set(time_step, displaced_surface);
1007
1008 // Eigen::MatrixXd forces = Eigen::MatrixXd::Zero(collision_mesh_.ndof(), 1);
1009 // forces -= barrier_potential_.gradient(collision_set, collision_mesh_, displaced_surface);
1010 // forces = collision_mesh_.to_full_dof(forces);
1011
1012 // double sum = (forces.array() * forces.array()).sum();
1013
1014 // return sum; },
1015 // G);
1016
1017 // logger().trace("gradu difference norm {}", (G - gradu).norm());
1018
1019 return weight() * gradu;
1020 }
1021
1022 void ProxyContactForceForm::compute_partial_gradient_step(const int time_step, const Eigen::VectorXd &x, Eigen::VectorXd &gradv) const
1023 {
1024 assert(state_.solve_data.time_integrator != nullptr);
1025 assert(state_.solve_data.contact_form != nullptr);
1026
1028 const Eigen::MatrixXd displaced_surface = collision_mesh_.displace_vertices(utils::unflatten(state_.diff_cached.u(time_step), collision_mesh_.dim()));
1029
1030 auto collision_set = get_or_compute_collision_set(time_step, displaced_surface);
1031
1032 Eigen::MatrixXd forces = collision_mesh_.to_full_dof(barrier_potential_.gradient(collision_set, collision_mesh_, displaced_surface));
1033
1034 StiffnessMatrix shape_derivative = collision_mesh_.to_full_dof(barrier_potential_.shape_derivative(collision_set, collision_mesh_, displaced_surface));
1035
1036 Eigen::VectorXd grads = 2 * shape_derivative.transpose() * forces;
1037
1038 // Eigen::VectorXd G;
1039 // fd::finite_gradient(
1040 // utils::flatten(node_positions_), [this, time_step](const Eigen::VectorXd &x_) {
1041 // Eigen::MatrixXd n;
1042 // Eigen::MatrixXi e, t;
1043 // Eigen::SparseMatrix<double> d_m;
1044 // std::vector<bool> i_o_s;
1045 // Eigen::MatrixXi c_c_c;
1046 // compute_forward_collision_mesh_quantities(state_, boundary_ids_, collision_mesh_, n, e, t, d_m, i_o_s, c_c_c);
1047
1048 // return compute_forward_quantity(
1049 // utils::unflatten(x_, n.cols()),
1050 // e,
1051 // t,
1052 // d_m,
1053 // i_o_s,
1054 // c_c_c,
1055 // dhat_,
1056 // dmin_,
1057 // [this, time_step](const Eigen::MatrixXd &displaced_surface) {return get_or_compute_collision_set(time_step, displaced_surface);},
1058 // state_.diff_cached.u(time_step),
1059 // barrier_potential_); },
1060 // G, fd::AccuracyOrder::SECOND, 1e-7);
1061
1062 // Eigen::MatrixXd diff(G.size(), 2);
1063 // diff.col(0) = G;
1064 // diff.col(1) = grads;
1065 // logger().trace("diff {}", diff);
1066 // logger().trace("size {} {}", G.size(), grads.size());
1067 // logger().trace("fd norm {}", G.norm());
1068 // logger().trace("grads difference norm {}", (G - grads).norm() / G.norm());
1069
1070 grads = state_.basis_nodes_to_gbasis_nodes * grads;
1071
1073 });
1074 }
1075
1076} // namespace polyfem::solver
ElementAssemblyValues vals
Definition Assembler.cpp:22
const double weight_
int x
main class that contains the polyfem solver and all its state
Definition State.hpp:79
int n_bases
number of bases
Definition State.hpp:178
int n_boundary_samples() const
quadrature used for projecting boundary conditions
Definition State.hpp:263
const std::vector< basis::ElementBases > & geom_bases() const
Get a constant reference to the geometry mapping bases.
Definition State.hpp:223
std::unique_ptr< mesh::Mesh > mesh
current mesh, it can be a Mesh2D or Mesh3D
Definition State.hpp:471
std::shared_ptr< assembler::Problem > problem
current problem, it contains rhs and bc
Definition State.hpp:168
json args
main input arguments containing all defaults
Definition State.hpp:101
solver::DiffCache diff_cached
Definition State.hpp:656
std::vector< basis::ElementBases > bases
FE bases, the size is #elements.
Definition State.hpp:171
std::vector< mesh::LocalBoundary > total_local_boundary
mapping from elements to nodes for all mesh
Definition State.hpp:436
solver::SolveData solve_data
timedependent stuff cached
Definition State.hpp:321
StiffnessMatrix basis_nodes_to_gbasis_nodes
Definition State.hpp:706
stores per local bases evaluations
std::vector< basis::Local2Global > global
stores per element basis values at given quadrature points and geometric mapping
void compute(const int el_index, const bool is_volume, const Eigen::MatrixXd &pts, const basis::ElementBases &basis, const basis::ElementBases &gbasis)
computes the per element values at the local (ref el) points (pts) sets basis_values,...
Stores the basis functions for a given element in a mesh (facet in 2d, cell in 3d).
Eigen::VectorXi local_nodes_for_primitive(const int local_index, const mesh::Mesh &mesh) const
static void extract_boundary_mesh(const mesh::Mesh &mesh, const int n_bases, const std::vector< basis::ElementBases > &bases, const std::vector< mesh::LocalBoundary > &total_local_boundary, Eigen::MatrixXd &node_positions, Eigen::MatrixXi &boundary_edges, Eigen::MatrixXi &boundary_triangles, std::vector< Eigen::Triplet< double > > &displacement_map_entries)
extracts the boundary mesh
Definition OutData.cpp:145
const VariableToSimulationGroup variable_to_simulations_
Eigen::VectorXd u(int step) const
virtual double weight() const
Get the form's multiplicative constant weight.
Definition Form.hpp:128
double value_unweighted_step(const int time_step, const Eigen::VectorXd &x) const override
std::vector< std::shared_ptr< ipc::NormalCollisions > > collision_sets_
std::map< int, std::set< int > > boundary_ids_to_dof_
Eigen::VectorXd compute_adjoint_rhs_step(const int time_step, const Eigen::VectorXd &x, const State &state) const override
void compute_partial_gradient_step(const int time_step, const Eigen::VectorXd &x, Eigen::VectorXd &gradv) const override
ProxyContactForceForm(const VariableToSimulationGroup &variable_to_simulations, const State &state, const double dhat, const bool quadratic_potential, const json &args)
const ipc::NormalCollisions & get_or_compute_collision_set(const int time_step, const Eigen::MatrixXd &displaced_surface) const
std::shared_ptr< solver::ContactForm > contact_form
std::shared_ptr< time_integrator::ImplicitTimeIntegrator > time_integrator
A collection of VariableToSimulation.
virtual Eigen::VectorXd apply_parametrization_jacobian(const ParameterType type, const State *state_ptr, const Eigen::VectorXd &x, const std::function< Eigen::VectorXd()> &grad) const
Maps the partial gradient wrt.
static bool boundary_quadrature(const mesh::LocalBoundary &local_boundary, const int order, const mesh::Mesh &mesh, const bool skip_computation, Eigen::MatrixXd &uv, Eigen::MatrixXd &points, Eigen::MatrixXd &normals, Eigen::VectorXd &weights, Eigen::VectorXi &global_primitive_ids)
str nodes
Definition p_bases.py:398
Eigen::VectorXd map_node_to_primitive_order(const State &state, const Eigen::VectorXd &nodes)
DScalar1< double, Eigen::Matrix< double, Eigen::Dynamic, 1 > > Diff
Eigen::Matrix< T, Eigen::Dynamic, Eigen::Dynamic, 0, 3, 3 > inverse(const Eigen::Matrix< T, Eigen::Dynamic, Eigen::Dynamic, 0, 3, 3 > &mat)
Eigen::MatrixXd unflatten(const Eigen::VectorXd &x, int dim)
Unflatten rowwises, so every dim elements in x become a row.
nlohmann::json json
Definition Common.hpp:9
void log_and_throw_adjoint_error(const std::string &msg)
Definition Logger.cpp:79
Eigen::SparseMatrix< double, Eigen::ColMajor > StiffnessMatrix
Definition Types.hpp:22
Automatic differentiation scalar with first-order derivatives.
Definition autodiff.h:112