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