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