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