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