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