PolyFEM
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
AdjointTools.cpp
Go to the documentation of this file.
1#include "AdjointTools.hpp"
2
3#include <polyfem/utils/MaybeParallelFor.hpp>
4#include <polyfem/utils/Timer.hpp>
5#include <polyfem/io/Evaluator.hpp>
6#include <polyfem/utils/AutodiffTypes.hpp>
7#include <polyfem/State.hpp>
8
9#include <polyfem/utils/IntegrableFunctional.hpp>
10#include <polyfem/utils/BoundarySampler.hpp>
11#include <polyfem/time_integrator/ImplicitTimeIntegrator.hpp>
12
13#include <polyfem/solver/forms/ElasticForm.hpp>
14#include <polyfem/solver/forms/ContactForm.hpp>
15#include <polyfem/solver/forms/BarrierContactForm.hpp>
16#include <polyfem/solver/forms/SmoothContactForm.hpp>
17#include <polyfem/solver/forms/PeriodicContactForm.hpp>
18#include <polyfem/solver/forms/NormalAdhesionForm.hpp>
19#include <polyfem/solver/forms/TangentialAdhesionForm.hpp>
20#include <polyfem/solver/forms/FrictionForm.hpp>
21#include <polyfem/solver/forms/BodyForm.hpp>
22#include <polyfem/solver/forms/PressureForm.hpp>
23#include <polyfem/solver/forms/InertiaForm.hpp>
24#include <polyfem/solver/forms/parametrization/PeriodicMeshToMesh.hpp>
25
26#include <polyfem/solver/NLProblem.hpp>
27#include <polyfem/solver/NLHomoProblem.hpp>
28
29#include <polyfem/time_integrator/BDF.hpp>
30
31/*
32Reminders:
33
34 1. Due to Dirichlet boundary, any force vector at dirichlet indices should be zero, so \partial_q h and \partial_u h should be set zero at dirichlet rows.
35
36*/
37
38using namespace polyfem::utils;
39
40namespace polyfem::solver
41{
42 namespace
43 {
44
45 int get_bdf_order(const polyfem::State &state)
46 {
47 if (state.args["time"]["integrator"].is_string())
48 return 1;
49 if (state.args["time"]["integrator"]["type"] == "ImplicitEuler")
50 return 1;
51 if (state.args["time"]["integrator"]["type"] == "BDF")
52 return state.args["time"]["integrator"]["steps"].get<int>();
53
54 polyfem::log_and_throw_adjoint_error("Integrator type not supported for differentiability.");
55 return -1;
56 }
57
58 double dot(const Eigen::MatrixXd &A, const Eigen::MatrixXd &B) { return (A.array() * B.array()).sum(); }
59
60 class LocalThreadScalarStorage
61 {
62 public:
63 double val;
64 assembler::ElementAssemblyValues vals;
65 QuadratureVector da;
66
67 LocalThreadScalarStorage()
68 {
69 val = 0;
70 }
71 };
72
73 class LocalThreadVecStorage
74 {
75 public:
76 Eigen::MatrixXd vec;
77 assembler::ElementAssemblyValues vals;
78 QuadratureVector da;
79
80 LocalThreadVecStorage(const int size)
81 {
82 vec.resize(size, 1);
83 vec.setZero();
84 }
85 };
86
87 typedef DScalar1<double, Eigen::Matrix<double, Eigen::Dynamic, 1>> Diff;
88
89 template <typename T>
90 T triangle_area(const Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> &V)
91 {
92 Eigen::Matrix<T, Eigen::Dynamic, 1> l1 = V.row(1) - V.row(0);
93 Eigen::Matrix<T, Eigen::Dynamic, 1> l2 = V.row(2) - V.row(0);
94 T area = 0.5 * sqrt(pow(l1(1) * l2(2) - l1(2) * l2(1), 2) + pow(l1(0) * l2(2) - l1(2) * l2(0), 2) + pow(l1(1) * l2(0) - l1(0) * l2(1), 2));
95 return area;
96 }
97
98 Eigen::MatrixXd triangle_area_grad(const Eigen::MatrixXd &F)
99 {
100 DiffScalarBase::setVariableCount(F.size());
101 Eigen::Matrix<Diff, Eigen::Dynamic, Eigen::Dynamic> full_diff(F.rows(), F.cols());
102 for (int i = 0; i < F.rows(); i++)
103 for (int j = 0; j < F.cols(); j++)
104 full_diff(i, j) = Diff(i + j * F.rows(), F(i, j));
105 auto reduced_diff = triangle_area(full_diff);
106
107 Eigen::MatrixXd grad(F.rows(), F.cols());
108 for (int i = 0; i < F.rows(); ++i)
109 for (int j = 0; j < F.cols(); ++j)
110 grad(i, j) = reduced_diff.getGradient()(i + j * F.rows());
111
112 return grad;
113 }
114
115 template <typename T>
116 T line_length(const Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> &V)
117 {
118 Eigen::Matrix<T, Eigen::Dynamic, 1> L = V.row(1) - V.row(0);
119 T area = L.norm();
120 return area;
121 }
122
123 Eigen::MatrixXd line_length_grad(const Eigen::MatrixXd &F)
124 {
125 DiffScalarBase::setVariableCount(F.size());
126 Eigen::Matrix<Diff, Eigen::Dynamic, Eigen::Dynamic> full_diff(F.rows(), F.cols());
127 for (int i = 0; i < F.rows(); i++)
128 for (int j = 0; j < F.cols(); j++)
129 full_diff(i, j) = Diff(i + j * F.rows(), F(i, j));
130 auto reduced_diff = line_length(full_diff);
131
132 Eigen::MatrixXd grad(F.rows(), F.cols());
133 for (int i = 0; i < F.rows(); ++i)
134 for (int j = 0; j < F.cols(); ++j)
135 grad(i, j) = reduced_diff.getGradient()(i + j * F.rows());
136
137 return grad;
138 }
139
140 template <typename T>
141 Eigen::Matrix<T, 2, 1> edge_normal(const Eigen::Matrix<T, 4, 1> &V)
142 {
143 Eigen::Matrix<T, 2, 1> v1 = V.segment(0, 2);
144 Eigen::Matrix<T, 2, 1> v2 = V.segment(2, 2);
145 Eigen::Matrix<T, 2, 1> normal = v1 - v2;
146 normal(0) *= -1;
147 normal = normal / normal.norm();
148 return normal;
149 }
150
151 template <typename T>
152 Eigen::Matrix<T, 3, 1> face_normal(const Eigen::Matrix<T, 9, 1> &V)
153 {
154 Eigen::Matrix<T, 3, 1> v1 = V.segment(0, 3);
155 Eigen::Matrix<T, 3, 1> v2 = V.segment(3, 3);
156 Eigen::Matrix<T, 3, 1> v3 = V.segment(6, 3);
157 Eigen::Matrix<T, 3, 1> normal = (v2 - v1).cross(v3 - v1);
158 normal = normal / normal.norm();
159 return normal;
160 }
161
162 Eigen::MatrixXd extract_lame_params(const std::map<std::string, Assembler::ParamFunc> &lame_params, const int e, const int t, const Eigen::MatrixXd &local_pts, const Eigen::MatrixXd &pts)
163 {
164 Eigen::MatrixXd params = Eigen::MatrixXd::Zero(local_pts.rows(), 2);
165
166 auto search_lambda = lame_params.find("lambda");
167 auto search_mu = lame_params.find("mu");
168
169 if (search_lambda == lame_params.end() || search_mu == lame_params.end())
170 return params;
171
172 for (int p = 0; p < local_pts.rows(); p++)
173 {
174 params(p, 0) = search_lambda->second(local_pts.row(p), pts.row(p), t, e);
175 params(p, 1) = search_mu->second(local_pts.row(p), pts.row(p), t, e);
176 }
177
178 return params;
179 }
180 } // namespace
181
182 double AdjointTools::integrate_objective(
183 const State &state,
184 const IntegrableFunctional &j,
185 const Eigen::MatrixXd &solution,
186 const std::set<int> &interested_ids, // either body id or surface id
187 const SpatialIntegralType spatial_integral_type,
188 const int cur_step) // current time step
189 {
190 const auto &bases = state.bases;
191 const auto &gbases = state.geom_bases();
192
193 const int dim = state.mesh->dimension();
194 const int actual_dim = state.problem->is_scalar() ? 1 : dim;
195 const int n_elements = int(bases.size());
196 const double t0 = state.problem->is_time_dependent() ? state.args["time"]["t0"].get<double>() : 0.0;
197 const double dt = state.problem->is_time_dependent() ? state.args["time"]["dt"].get<double>() : 0.0;
198
199 double integral = 0;
200 if (spatial_integral_type == SpatialIntegralType::Volume)
201 {
202 auto storage = utils::create_thread_storage(LocalThreadScalarStorage());
203 utils::maybe_parallel_for(n_elements, [&](int start, int end, int thread_id) {
204 LocalThreadScalarStorage &local_storage = utils::get_local_thread_storage(storage, thread_id);
205
206 IntegrableFunctional::ParameterType params;
207 params.t = dt * cur_step + t0;
208 params.step = cur_step;
209
210 Eigen::MatrixXd u, grad_u;
211 Eigen::MatrixXd result;
212
213 for (int e = start; e < end; ++e)
214 {
215 if (interested_ids.size() != 0 && interested_ids.find(state.mesh->get_body_id(e)) == interested_ids.end())
216 continue;
217
218 assembler::ElementAssemblyValues &vals = local_storage.vals;
219 state.ass_vals_cache.compute(e, state.mesh->is_volume(), bases[e], gbases[e], vals);
220 io::Evaluator::interpolate_at_local_vals(e, dim, actual_dim, vals, solution, u, grad_u);
221
222 const quadrature::Quadrature &quadrature = vals.quadrature;
223 local_storage.da = vals.det.array() * quadrature.weights.array();
224
225 const Eigen::MatrixXd lame_params = extract_lame_params(state.assembler->parameters(), e, params.t, quadrature.points, vals.val);
226
227 params.elem = e;
228 params.body_id = state.mesh->get_body_id(e);
229 j.evaluate(lame_params, quadrature.points, vals.val, u, grad_u, Eigen::MatrixXd::Zero(0, 0) /*Not used*/, vals, params, result);
230
231 local_storage.val += dot(result, local_storage.da);
232 }
233 });
234 for (const LocalThreadScalarStorage &local_storage : storage)
235 integral += local_storage.val;
236 }
237 else if (spatial_integral_type == SpatialIntegralType::Surface)
238 {
239 auto storage = utils::create_thread_storage(LocalThreadScalarStorage());
240 utils::maybe_parallel_for(state.total_local_boundary.size(), [&](int start, int end, int thread_id) {
241 LocalThreadScalarStorage &local_storage = utils::get_local_thread_storage(storage, thread_id);
242
243 Eigen::MatrixXd uv;
244 Eigen::MatrixXd points, normal;
245 Eigen::VectorXd weights;
246
247 Eigen::MatrixXd u, grad_u;
248 Eigen::MatrixXd result;
249 IntegrableFunctional::ParameterType params;
250 params.t = dt * cur_step + t0;
251 params.step = cur_step;
252
253 for (int lb_id = start; lb_id < end; ++lb_id)
254 {
255 const auto &lb = state.total_local_boundary[lb_id];
256 const int e = lb.element_id();
257
258 for (int i = 0; i < lb.size(); i++)
259 {
260 const int global_primitive_id = lb.global_primitive_id(i);
261 if (interested_ids.size() != 0 && interested_ids.find(state.mesh->get_boundary_id(global_primitive_id)) == interested_ids.end())
262 continue;
263
264 utils::BoundarySampler::boundary_quadrature(lb, state.n_boundary_samples(), *state.mesh, i, false, uv, points, normal, weights);
265
266 assembler::ElementAssemblyValues &vals = local_storage.vals;
267 vals.compute(e, state.mesh->is_volume(), points, bases[e], gbases[e]);
268 io::Evaluator::interpolate_at_local_vals(e, dim, actual_dim, vals, solution, u, grad_u);
269
270 const Eigen::MatrixXd lame_params = extract_lame_params(state.assembler->parameters(), e, params.t, points, vals.val);
271
272 params.elem = e;
273 params.body_id = state.mesh->get_body_id(e);
274 params.boundary_id = state.mesh->get_boundary_id(global_primitive_id);
275 j.evaluate(lame_params, points, vals.val, u, grad_u, normal, vals, params, result);
276
277 local_storage.val += dot(result, weights);
278 }
279 }
280 });
281 for (const LocalThreadScalarStorage &local_storage : storage)
282 integral += local_storage.val;
283 }
284 else if (spatial_integral_type == SpatialIntegralType::VertexSum)
285 {
286 std::vector<bool> traversed(state.n_bases, false);
287 IntegrableFunctional::ParameterType params;
288 params.t = dt * cur_step + t0;
289 params.step = cur_step;
290 for (int e = 0; e < bases.size(); e++)
291 {
292 const auto &bs = bases[e];
293 for (int i = 0; i < bs.bases.size(); i++)
294 {
295 const auto &b = bs.bases[i];
296 assert(b.global().size() == 1);
297 const auto &g = b.global()[0];
298 if (traversed[g.index])
299 continue;
300
301 const Eigen::MatrixXd lame_params = extract_lame_params(state.assembler->parameters(), e, params.t, Eigen::MatrixXd::Zero(1, dim) /*Not used*/, g.node);
302
303 params.node = g.index;
304 params.elem = e;
305 params.body_id = state.mesh->get_body_id(e);
306 Eigen::MatrixXd val;
307 j.evaluate(lame_params, Eigen::MatrixXd::Zero(1, dim) /*Not used*/, g.node, solution.block(g.index * dim, 0, dim, 1).transpose(), Eigen::MatrixXd::Zero(1, dim * actual_dim) /*Not used*/, Eigen::MatrixXd::Zero(0, 0) /*Not used*/, assembler::ElementAssemblyValues(), params, val);
308 integral += val(0);
309 traversed[g.index] = true;
310 }
311 }
312 }
313
314 return integral;
315 }
316
317 void AdjointTools::compute_shape_derivative_functional_term(
318 const State &state,
319 const Eigen::MatrixXd &solution,
320 const IntegrableFunctional &j,
321 const std::set<int> &interested_ids, // either body id or surface id
322 const SpatialIntegralType spatial_integral_type,
323 Eigen::VectorXd &term,
324 const int cur_time_step)
325 {
326 const auto &gbases = state.geom_bases();
327 const auto &bases = state.bases;
328 const int dim = state.mesh->dimension();
329 const int actual_dim = state.problem->is_scalar() ? 1 : dim;
330 const double t0 = state.problem->is_time_dependent() ? state.args["time"]["t0"].get<double>() : 0.0;
331 const double dt = state.problem->is_time_dependent() ? state.args["time"]["dt"].get<double>() : 0.0;
332
333 const int n_elements = int(bases.size());
334 term.setZero(state.n_geom_bases * dim, 1);
335
336 auto storage = utils::create_thread_storage(LocalThreadVecStorage(term.size()));
337
338 if (spatial_integral_type == SpatialIntegralType::Volume)
339 {
340 utils::maybe_parallel_for(n_elements, [&](int start, int end, int thread_id) {
341 LocalThreadVecStorage &local_storage = utils::get_local_thread_storage(storage, thread_id);
342
343 Eigen::MatrixXd u, grad_u, j_val, dj_dgradu, dj_dx;
344
345 IntegrableFunctional::ParameterType params;
346 params.t = cur_time_step * dt + t0;
347 params.step = cur_time_step;
348
349 for (int e = start; e < end; ++e)
350 {
351 if (interested_ids.size() != 0 && interested_ids.find(state.mesh->get_body_id(e)) == interested_ids.end())
352 continue;
353
354 assembler::ElementAssemblyValues &vals = local_storage.vals;
355 state.ass_vals_cache.compute(e, state.mesh->is_volume(), bases[e], gbases[e], vals);
356 io::Evaluator::interpolate_at_local_vals(e, dim, actual_dim, vals, solution, u, grad_u);
357
358 assembler::ElementAssemblyValues gvals;
359 gvals.compute(e, state.mesh->is_volume(), vals.quadrature.points, gbases[e], gbases[e]);
360
361 const quadrature::Quadrature &quadrature = vals.quadrature;
362 local_storage.da = vals.det.array() * quadrature.weights.array();
363
364 const Eigen::MatrixXd lame_params = extract_lame_params(state.assembler->parameters(), e, params.t, quadrature.points, vals.val);
365
366 params.elem = e;
367 params.body_id = state.mesh->get_body_id(e);
368
369 j.evaluate(lame_params, quadrature.points, vals.val, u, grad_u, Eigen::MatrixXd::Zero(0, 0) /*Not used*/, vals, params, j_val);
370
371 if (j.depend_on_gradu())
372 j.dj_dgradu(lame_params, quadrature.points, vals.val, u, grad_u, Eigen::MatrixXd::Zero(0, 0) /*Not used*/, vals, params, dj_dgradu);
373
374 if (j.depend_on_x())
375 j.dj_dx(lame_params, quadrature.points, vals.val, u, grad_u, Eigen::MatrixXd::Zero(0, 0) /*Not used*/, vals, params, dj_dx);
376
377 Eigen::MatrixXd tau_q, grad_u_q;
378 for (auto &v : gvals.basis_values)
379 {
380 for (int q = 0; q < local_storage.da.size(); ++q)
381 {
382 local_storage.vec.block(v.global[0].index * dim, 0, dim, 1) += (j_val(q) * local_storage.da(q)) * v.grad_t_m.row(q).transpose();
383
384 if (j.depend_on_x())
385 local_storage.vec.block(v.global[0].index * dim, 0, dim, 1) += (v.val(q) * local_storage.da(q)) * dj_dx.row(q).transpose();
386
387 if (j.depend_on_gradu())
388 {
389 if (dim == actual_dim) // Elasticity PDE
390 {
391 vector2matrix(dj_dgradu.row(q), tau_q);
392 vector2matrix(grad_u.row(q), grad_u_q);
393 }
394 else // Laplacian PDE
395 {
396 tau_q = dj_dgradu.row(q);
397 grad_u_q = grad_u.row(q);
398 }
399 for (int d = 0; d < dim; d++)
400 local_storage.vec(v.global[0].index * dim + d) += -dot(tau_q, grad_u_q.col(d) * v.grad_t_m.row(q)) * local_storage.da(q);
401 }
402 }
403 }
404 }
405 });
406 }
407 else if (spatial_integral_type == SpatialIntegralType::Surface)
408 {
409 utils::maybe_parallel_for(state.total_local_boundary.size(), [&](int start, int end, int thread_id) {
410 LocalThreadVecStorage &local_storage = utils::get_local_thread_storage(storage, thread_id);
411
412 Eigen::MatrixXd uv, points, normal;
413 Eigen::VectorXd &weights = local_storage.da;
414
415 Eigen::MatrixXd u, grad_u, x, grad_x, j_val, dj_dgradu, dj_dgradx, dj_dx;
416
417 IntegrableFunctional::ParameterType params;
418 params.t = cur_time_step * dt + t0;
419 params.step = cur_time_step;
420
421 for (int lb_id = start; lb_id < end; ++lb_id)
422 {
423 const auto &lb = state.total_local_boundary[lb_id];
424 const int e = lb.element_id();
425
426 for (int i = 0; i < lb.size(); i++)
427 {
428 const int global_primitive_id = lb.global_primitive_id(i);
429 if (interested_ids.size() != 0 && interested_ids.find(state.mesh->get_boundary_id(global_primitive_id)) == interested_ids.end())
430 continue;
431
432 utils::BoundarySampler::boundary_quadrature(lb, state.n_boundary_samples(), *state.mesh, i, false, uv, points, normal, weights);
433
434 assembler::ElementAssemblyValues &vals = local_storage.vals;
435 io::Evaluator::interpolate_at_local_vals(*state.mesh, state.problem->is_scalar(), bases, gbases, e, points, solution, u, grad_u);
436 // io::Evaluator::interpolate_at_local_vals(*state.mesh, state.problem->is_scalar(), gbases, gbases, e, points, global_positions, x, grad_x);
437
438 vals.compute(e, state.mesh->is_volume(), points, gbases[e], gbases[e]);
439
440 // normal = normal * vals.jac_it[0]; // assuming linear geometry
441
442 const int n_loc_bases_ = int(vals.basis_values.size());
443
444 const Eigen::MatrixXd lame_params = extract_lame_params(state.assembler->parameters(), e, params.t, points, vals.val);
445
446 params.elem = e;
447 params.body_id = state.mesh->get_body_id(e);
448 params.boundary_id = state.mesh->get_boundary_id(global_primitive_id);
449
450 j.evaluate(lame_params, points, vals.val, u, grad_u, normal, vals, params, j_val);
451 j_val = j_val.array().colwise() * weights.array();
452
453 if (j.depend_on_gradu())
454 {
455 j.dj_dgradu(lame_params, points, vals.val, u, grad_u, normal, vals, params, dj_dgradu);
456 dj_dgradu = dj_dgradu.array().colwise() * weights.array();
457 }
458
459 if (j.depend_on_gradx())
460 {
461 j.dj_dgradx(lame_params, points, vals.val, u, grad_u, normal, vals, params, dj_dgradx);
462 dj_dgradx = dj_dgradx.array().colwise() * weights.array();
463 }
464
465 if (j.depend_on_x())
466 {
467 j.dj_dx(lame_params, points, vals.val, u, grad_u, normal, vals, params, dj_dx);
468 dj_dx = dj_dx.array().colwise() * weights.array();
469 }
470
471 const auto nodes = gbases[e].local_nodes_for_primitive(lb.global_primitive_id(i), *state.mesh);
472
473 if (nodes.size() != dim)
474 log_and_throw_adjoint_error("Only linear geometry is supported in differentiable surface integral functional!");
475
476 Eigen::MatrixXd velocity_div_mat;
477 if (state.mesh->is_volume())
478 {
479 Eigen::Matrix3d V;
480 for (int d = 0; d < 3; d++)
481 V.row(d) = gbases[e].bases[nodes(d)].global()[0].node;
482 velocity_div_mat = face_velocity_divergence(V);
483 }
484 else
485 {
486 Eigen::Matrix2d V;
487 for (int d = 0; d < 2; d++)
488 V.row(d) = gbases[e].bases[nodes(d)].global()[0].node;
489 velocity_div_mat = edge_velocity_divergence(V);
490 }
491
492 Eigen::MatrixXd grad_u_q, tau_q, grad_x_q;
493 for (long n = 0; n < nodes.size(); ++n)
494 {
495 const assembler::AssemblyValues &v = vals.basis_values[nodes(n)];
496
497 local_storage.vec.block(v.global[0].index * dim, 0, dim, 1) += j_val.sum() * velocity_div_mat.row(n).transpose();
498 }
499
500 for (long n = 0; n < n_loc_bases_; ++n)
501 {
502 const assembler::AssemblyValues &v = vals.basis_values[n];
503
504 if (j.depend_on_x())
505 local_storage.vec.block(v.global[0].index * dim, 0, dim, 1) += dj_dx.transpose() * v.val;
506
507 // integrate j * div(gbases) over the whole boundary
508 if (j.depend_on_gradu())
509 {
510 for (int q = 0; q < weights.size(); ++q)
511 {
512 if (dim == actual_dim) // Elasticity PDE
513 {
514 vector2matrix(grad_u.row(q), grad_u_q);
515 vector2matrix(dj_dgradu.row(q), tau_q);
516 }
517 else // Laplacian PDE
518 {
519 grad_u_q = grad_u.row(q);
520 tau_q = dj_dgradu.row(q);
521 }
522
523 for (int d = 0; d < dim; d++)
524 local_storage.vec(v.global[0].index * dim + d) += -dot(tau_q, grad_u_q.col(d) * v.grad_t_m.row(q));
525 }
526 }
527
528 if (j.depend_on_gradx())
529 {
530 for (int d = 0; d < dim; d++)
531 {
532 for (int q = 0; q < weights.size(); ++q)
533 local_storage.vec(v.global[0].index * dim + d) += dot(dj_dgradx.block(q, d * dim, 1, dim), v.grad.row(q));
534 }
535 }
536 }
537 }
538 }
539 });
540 }
541 else if (spatial_integral_type == SpatialIntegralType::VertexSum)
542 {
543 log_and_throw_adjoint_error("Shape derivative of vertex sum type functional is not implemented!");
544 }
545 for (const LocalThreadVecStorage &local_storage : storage)
546 term += local_storage.vec;
547
548 term = utils::flatten(utils::unflatten(term, dim)(state.primitive_to_node(), Eigen::all));
549 }
550
551 void AdjointTools::dJ_shape_static_adjoint_term(
552 const State &state,
553 const Eigen::MatrixXd &sol,
554 const Eigen::MatrixXd &adjoint,
555 Eigen::VectorXd &one_form)
556 {
557 Eigen::VectorXd elasticity_term, rhs_term, pressure_term, contact_term, adhesion_term;
558
559 one_form.setZero(state.n_geom_bases * state.mesh->dimension());
560 Eigen::MatrixXd adjoint_zeroed = adjoint;
561 adjoint_zeroed(state.boundary_nodes, Eigen::all).setZero();
562
563 // if (j.depend_on_u() || j.depend_on_gradu())
564 {
565 state.solve_data.elastic_form->force_shape_derivative(0, state.n_geom_bases, sol, sol, adjoint_zeroed, elasticity_term);
566 if (state.solve_data.body_form)
567 state.solve_data.body_form->force_shape_derivative(state.n_geom_bases, 0, sol, adjoint_zeroed, rhs_term);
568 else
569 rhs_term.setZero(one_form.size());
570
571 if (state.solve_data.pressure_form)
572 {
573 state.solve_data.pressure_form->force_shape_derivative(state.n_geom_bases, 0, sol, adjoint_zeroed, pressure_term);
574 pressure_term = state.basis_nodes_to_gbasis_nodes * pressure_term;
575 }
576 else
577 pressure_term.setZero(one_form.size());
578
579 if (state.is_contact_enabled())
580 {
581 if (const auto barrier_contact = dynamic_cast<const BarrierContactForm*>(state.solve_data.contact_form.get()))
582 {
583 barrier_contact->force_shape_derivative(state.diff_cached.collision_set(0), sol, adjoint_zeroed, contact_term);
584 }
585 else if (const auto smooth_contact = dynamic_cast<const SmoothContactForm*>(state.solve_data.contact_form.get()))
586 {
587 smooth_contact->force_shape_derivative(state.diff_cached.smooth_collision_set(0), sol, adjoint_zeroed, contact_term);
588 }
589
590 contact_term = state.basis_nodes_to_gbasis_nodes * contact_term;
591 }
592 else
593 contact_term.setZero(elasticity_term.size());
594
595 if (state.is_adhesion_enabled())
596 {
597 state.solve_data.normal_adhesion_form->force_shape_derivative(state.diff_cached.normal_adhesion_collision_set(0), sol, adjoint, adhesion_term);
598 adhesion_term = state.basis_nodes_to_gbasis_nodes * adhesion_term;
599 }
600 else
601 {
602 adhesion_term.setZero(elasticity_term.size());
603 }
604 }
605
606 one_form -= elasticity_term + rhs_term + pressure_term + contact_term + adhesion_term;
607 one_form = utils::flatten(utils::unflatten(one_form, state.mesh->dimension())(state.primitive_to_node(), Eigen::all));
608 }
609
610 void AdjointTools::dJ_shape_homogenization_adjoint_term(
611 const State &state,
612 const Eigen::MatrixXd &sol,
613 const Eigen::MatrixXd &adjoint,
614 Eigen::VectorXd &one_form)
615 {
616 Eigen::VectorXd elasticity_term, contact_term, adhesion_term;
617
618 std::shared_ptr<NLHomoProblem> homo_problem = std::dynamic_pointer_cast<NLHomoProblem>(state.solve_data.nl_problem);
619 assert(homo_problem);
620
621 const int dim = state.mesh->dimension();
622 one_form.setZero(state.n_geom_bases * dim);
623
624 const Eigen::MatrixXd affine_adjoint = homo_problem->reduced_to_disp_grad(adjoint, true);
625 const Eigen::VectorXd full_adjoint = homo_problem->NLProblem::reduced_to_full(adjoint.topRows(homo_problem->reduced_size())) + io::Evaluator::generate_linear_field(state.n_bases, state.mesh_nodes, affine_adjoint);
626
627 state.solve_data.elastic_form->force_shape_derivative(0, state.n_geom_bases, sol, sol, full_adjoint, elasticity_term);
628
629 if (state.solve_data.contact_form)
630 {
631 if (const auto barrier_contact = dynamic_cast<const BarrierContactForm*>(state.solve_data.contact_form.get()))
632 {
633 barrier_contact->force_shape_derivative(state.diff_cached.collision_set(0), sol, full_adjoint, contact_term);
634 }
635 else if (const auto smooth_contact = dynamic_cast<const SmoothContactForm*>(state.solve_data.contact_form.get()))
636 {
637 smooth_contact->force_shape_derivative(state.diff_cached.smooth_collision_set(0), sol, full_adjoint, contact_term);
638 }
639
640 contact_term = state.basis_nodes_to_gbasis_nodes * contact_term;
641 }
642 else
643 contact_term.setZero(elasticity_term.size());
644
645 if (state.is_adhesion_enabled())
646 {
647 state.solve_data.normal_adhesion_form->force_shape_derivative(state.diff_cached.normal_adhesion_collision_set(0), sol, full_adjoint, adhesion_term);
648 adhesion_term = state.basis_nodes_to_gbasis_nodes * adhesion_term;
649 }
650 else
651 {
652 adhesion_term.setZero(elasticity_term.size());
653 }
654
655 one_form = -(elasticity_term + contact_term + adhesion_term);
656
657 Eigen::VectorXd force;
658 homo_problem->FullNLProblem::gradient(sol, force);
659 one_form -= state.basis_nodes_to_gbasis_nodes * utils::flatten(utils::unflatten(force, dim) * affine_adjoint);
660
661 one_form = utils::flatten(utils::unflatten(one_form, dim)(state.primitive_to_node(), Eigen::all));
662 }
663
664 void AdjointTools::dJ_periodic_shape_adjoint_term(
665 const State &state,
666 const PeriodicMeshToMesh &periodic_mesh_map,
667 const Eigen::VectorXd &periodic_mesh_representation,
668 const Eigen::MatrixXd &sol,
669 const Eigen::MatrixXd &adjoint,
670 Eigen::VectorXd &one_form)
671 {
672 std::shared_ptr<NLHomoProblem> homo_problem = std::dynamic_pointer_cast<NLHomoProblem>(state.solve_data.nl_problem);
673 assert(homo_problem);
674
675 const Eigen::MatrixXd reduced_sol = homo_problem->full_to_reduced(sol, state.diff_cached.disp_grad());
676 const Eigen::VectorXd extended_sol = homo_problem->reduced_to_extended(reduced_sol);
677
678 const Eigen::VectorXd extended_adjoint = homo_problem->reduced_to_extended(adjoint, true);
679 const Eigen::MatrixXd affine_adjoint = homo_problem->reduced_to_disp_grad(adjoint, true);
680 const Eigen::VectorXd full_adjoint = homo_problem->NLProblem::reduced_to_full(adjoint.topRows(homo_problem->reduced_size())) + io::Evaluator::generate_linear_field(state.n_bases, state.mesh_nodes, affine_adjoint);
681
682 const int dim = state.mesh->dimension();
683
684 dJ_shape_homogenization_adjoint_term(state, sol, adjoint, one_form);
685
686 StiffnessMatrix hessian;
687 homo_problem->set_project_to_psd(false);
688 homo_problem->FullNLProblem::hessian(sol, hessian);
689 Eigen::VectorXd partial_term = full_adjoint.transpose() * hessian;
690 partial_term = state.basis_nodes_to_gbasis_nodes * utils::flatten(utils::unflatten(partial_term, dim) * state.diff_cached.disp_grad());
691 one_form -= utils::flatten(utils::unflatten(partial_term, dim)(state.primitive_to_node(), Eigen::all));
692
693 one_form = periodic_mesh_map.apply_jacobian(one_form, periodic_mesh_representation);
694
695 if (state.solve_data.periodic_contact_form)
696 {
697 Eigen::VectorXd contact_term;
698 state.solve_data.periodic_contact_form->force_periodic_shape_derivative(state, periodic_mesh_map, periodic_mesh_representation, state.solve_data.periodic_contact_form->collision_set(), extended_sol, extended_adjoint, contact_term);
699
700 one_form -= contact_term;
701 }
702 }
703
704 void AdjointTools::dJ_shape_transient_adjoint_term(
705 const State &state,
706 const Eigen::MatrixXd &adjoint_nu,
707 const Eigen::MatrixXd &adjoint_p,
708 Eigen::VectorXd &one_form)
709 {
710 const double t0 = state.args["time"]["t0"];
711 const double dt = state.args["time"]["dt"];
712 const int time_steps = state.args["time"]["time_steps"];
713 const int bdf_order = get_bdf_order(state);
714
715 Eigen::VectorXd elasticity_term, rhs_term, pressure_term, damping_term, mass_term, contact_term, friction_term, adhesion_term, tangential_adhesion_term;
716 one_form.setZero(state.n_geom_bases * state.mesh->dimension());
717
718 Eigen::VectorXd cur_p, cur_nu;
719 for (int i = time_steps; i > 0; --i)
720 {
721 const int real_order = std::min(bdf_order, i);
722 double beta = time_integrator::BDF::betas(real_order - 1);
723 double beta_dt = beta * dt;
724 const double t = i * dt + t0;
725
726 Eigen::MatrixXd velocity = state.diff_cached.v(i);
727
728 cur_p = adjoint_p.col(i);
729 cur_nu = adjoint_nu.col(i);
730 cur_p(state.boundary_nodes).setZero();
731 cur_nu(state.boundary_nodes).setZero();
732
733 {
734 state.solve_data.inertia_form->force_shape_derivative(state.mesh->is_volume(), state.n_geom_bases, t, state.bases, state.geom_bases(), *(state.mass_matrix_assembler), state.mass_ass_vals_cache, velocity, cur_nu, mass_term);
735 state.solve_data.elastic_form->force_shape_derivative(t, state.n_geom_bases, state.diff_cached.u(i), state.diff_cached.u(i), cur_p, elasticity_term);
736 state.solve_data.body_form->force_shape_derivative(state.n_geom_bases, t, state.diff_cached.u(i - 1), cur_p, rhs_term);
737 state.solve_data.pressure_form->force_shape_derivative(state.n_geom_bases, t, state.diff_cached.u(i), cur_p, pressure_term);
738 pressure_term = state.basis_nodes_to_gbasis_nodes * pressure_term;
739
740 if (state.solve_data.damping_form)
741 state.solve_data.damping_form->force_shape_derivative(t, state.n_geom_bases, state.diff_cached.u(i), state.diff_cached.u(i - 1), cur_p, damping_term);
742 else
743 damping_term.setZero(mass_term.size());
744
745 if (state.is_contact_enabled())
746 {
747 if (const auto barrier_contact = dynamic_cast<const BarrierContactForm*>(state.solve_data.contact_form.get()))
748 {
749 barrier_contact->force_shape_derivative(state.diff_cached.collision_set(i), state.diff_cached.u(i), cur_p, contact_term);
750 }
751 else if (const auto smooth_contact = dynamic_cast<const SmoothContactForm*>(state.solve_data.contact_form.get()))
752 {
753 smooth_contact->force_shape_derivative(state.diff_cached.smooth_collision_set(i), state.diff_cached.u(i), cur_p, contact_term);
754 }
755 contact_term = state.basis_nodes_to_gbasis_nodes * contact_term;
756 // contact_term /= beta_dt * beta_dt;
757 }
758 else
759 contact_term.setZero(mass_term.size());
760
761 if (state.solve_data.friction_form)
762 {
763 state.solve_data.friction_form->force_shape_derivative(state.diff_cached.u(i - 1), state.diff_cached.u(i), cur_p, state.diff_cached.friction_collision_set(i), friction_term);
764 friction_term = state.basis_nodes_to_gbasis_nodes * (friction_term / beta);
765 // friction_term /= beta_dt * beta_dt;
766 }
767 else
768 friction_term.setZero(mass_term.size());
769
770 if (state.is_adhesion_enabled())
771 {
772 state.solve_data.normal_adhesion_form->force_shape_derivative(state.diff_cached.normal_adhesion_collision_set(i), state.diff_cached.u(i), cur_p, adhesion_term);
773 adhesion_term = state.basis_nodes_to_gbasis_nodes * adhesion_term;
774 }
775 else
776 {
777 adhesion_term.setZero(mass_term.size());
778 }
779
780 if (state.solve_data.tangential_adhesion_form)
781 {
782 state.solve_data.tangential_adhesion_form->force_shape_derivative(state.diff_cached.u(i - 1), state.diff_cached.u(i), cur_p, state.diff_cached.tangential_adhesion_collision_set(i), tangential_adhesion_term);
783 tangential_adhesion_term = state.basis_nodes_to_gbasis_nodes * (tangential_adhesion_term / beta);
784 // friction_term /= beta_dt * beta_dt;
785 }
786 else
787 tangential_adhesion_term.setZero(mass_term.size());
788 }
789
790 one_form += beta_dt * (elasticity_term + rhs_term + pressure_term + damping_term + contact_term + friction_term + mass_term + adhesion_term + tangential_adhesion_term);
791 }
792
793 // time step 0
794 Eigen::VectorXd sum_alpha_p;
795 {
796 sum_alpha_p.setZero(adjoint_p.rows());
797 int num = std::min(bdf_order, time_steps);
798 for (int j = 0; j < num; ++j)
799 {
800 int order = std::min(bdf_order - 1, j);
801 sum_alpha_p -= time_integrator::BDF::alphas(order)[j] * adjoint_p.col(j + 1);
802 }
803 }
804 sum_alpha_p(state.boundary_nodes).setZero();
805 state.solve_data.inertia_form->force_shape_derivative(state.mesh->is_volume(), state.n_geom_bases, t0, state.bases, state.geom_bases(), *(state.mass_matrix_assembler), state.mass_ass_vals_cache, state.diff_cached.v(0), sum_alpha_p, mass_term);
806
807 one_form += mass_term;
808
809 one_form = utils::flatten(utils::unflatten(one_form, state.mesh->dimension())(state.primitive_to_node(), Eigen::all));
810 }
811
812 void AdjointTools::dJ_material_static_adjoint_term(
813 const State &state,
814 const Eigen::MatrixXd &sol,
815 const Eigen::MatrixXd &adjoint,
816 Eigen::VectorXd &one_form)
817 {
818 Eigen::MatrixXd adjoint_zeroed = adjoint;
819 adjoint_zeroed(state.boundary_nodes, Eigen::all).setZero();
820 state.solve_data.elastic_form->force_material_derivative(0, sol, sol, adjoint_zeroed, one_form);
821 }
822
823 void AdjointTools::dJ_material_transient_adjoint_term(
824 const State &state,
825 const Eigen::MatrixXd &adjoint_nu,
826 const Eigen::MatrixXd &adjoint_p,
827 Eigen::VectorXd &one_form)
828 {
829 const double t0 = state.args["time"]["t0"];
830 const double dt = state.args["time"]["dt"];
831 const int time_steps = state.args["time"]["time_steps"];
832 const int bdf_order = get_bdf_order(state);
833
834 one_form.setZero(state.bases.size() * 2);
835
836 auto storage = utils::create_thread_storage(LocalThreadVecStorage(one_form.size()));
837
838 utils::maybe_parallel_for(time_steps, [&](int start, int end, int thread_id) {
839 LocalThreadVecStorage &local_storage = utils::get_local_thread_storage(storage, thread_id);
840 Eigen::VectorXd elasticity_term;
841 for (int i_aux = start; i_aux < end; ++i_aux)
842 {
843 const int i = time_steps - i_aux;
844 const int real_order = std::min(bdf_order, i);
845 double beta_dt = time_integrator::BDF::betas(real_order - 1) * dt;
846
847 Eigen::VectorXd cur_p = adjoint_p.col(i);
848 cur_p(state.boundary_nodes).setZero();
849
850 state.solve_data.elastic_form->force_material_derivative(t0 + dt * i, state.diff_cached.u(i), state.diff_cached.u(i - 1), -cur_p, elasticity_term);
851 local_storage.vec += beta_dt * elasticity_term;
852 }
853 });
854
855 for (const LocalThreadVecStorage &local_storage : storage)
856 one_form += local_storage.vec;
857 }
858
859 void AdjointTools::dJ_friction_transient_adjoint_term(
860 const State &state,
861 const Eigen::MatrixXd &adjoint_nu,
862 const Eigen::MatrixXd &adjoint_p,
863 Eigen::VectorXd &one_form)
864 {
865 const double dt = state.args["time"]["dt"];
866 const double mu = state.solve_data.friction_form->mu();
867 const int time_steps = state.args["time"]["time_steps"];
868 const int dim = state.mesh->dimension();
869 const int bdf_order = get_bdf_order(state);
870
871 one_form.setZero(1);
872
873 std::shared_ptr<time_integrator::ImplicitTimeIntegrator> time_integrator =
874 time_integrator::ImplicitTimeIntegrator::construct_time_integrator(state.args["time"]["integrator"]);
875 {
876 Eigen::MatrixXd solution, velocity, acceleration;
877 solution = state.diff_cached.u(0);
878 state.initial_velocity(velocity);
879 state.initial_acceleration(acceleration);
880 if (state.initial_vel_update.size() == state.ndof())
881 velocity = state.initial_vel_update;
882 const double dt = state.args["time"]["dt"];
883 time_integrator->init(solution, velocity, acceleration, dt);
884 }
885
886 for (int t = 1; t <= time_steps; ++t)
887 {
888 const int real_order = std::min(bdf_order, t);
889 double beta = time_integrator::BDF::betas(real_order - 1);
890
891 const Eigen::MatrixXd surface_solution_prev = state.collision_mesh.vertices(utils::unflatten(state.diff_cached.u(t - 1), dim));
892 // const Eigen::MatrixXd surface_solution = state.collision_mesh.vertices(utils::unflatten(state.diff_cached.u(t), dim));
893
894 const Eigen::MatrixXd surface_velocities = state.collision_mesh.map_displacements(utils::unflatten(time_integrator->compute_velocity(state.diff_cached.u(t)), state.collision_mesh.dim()));
895 time_integrator->update_quantities(state.diff_cached.u(t));
896
897 if (const auto barrier_contact = dynamic_cast<const BarrierContactForm*>(state.solve_data.contact_form.get()))
898 {
899 Eigen::MatrixXd force = state.collision_mesh.to_full_dof(
900 -state.solve_data.friction_form->friction_potential().force(
901 state.diff_cached.friction_collision_set(t),
902 state.collision_mesh,
903 state.collision_mesh.rest_positions(),
904 /*lagged_displacements=*/surface_solution_prev,
905 surface_velocities,
906 barrier_contact->barrier_potential(),
907 barrier_contact->barrier_stiffness(),
908 0., true));
909
910 Eigen::VectorXd cur_p = adjoint_p.col(t);
911 cur_p(state.boundary_nodes).setZero();
912
913 one_form(0) += dot(cur_p, force) * beta * dt;
914 }
915 }
916 }
917
918 void AdjointTools::dJ_damping_transient_adjoint_term(
919 const State &state,
920 const Eigen::MatrixXd &adjoint_nu,
921 const Eigen::MatrixXd &adjoint_p,
922 Eigen::VectorXd &one_form)
923 {
924 const double t0 = state.args["time"]["t0"];
925 const double dt = state.args["time"]["dt"];
926 const int time_steps = state.args["time"]["time_steps"];
927 const int bdf_order = get_bdf_order(state);
928
929 one_form.setZero(2);
930
931 auto storage = utils::create_thread_storage(LocalThreadVecStorage(one_form.size()));
932
933 utils::maybe_parallel_for(time_steps, [&](int start, int end, int thread_id) {
934 LocalThreadVecStorage &local_storage = utils::get_local_thread_storage(storage, thread_id);
935 Eigen::VectorXd damping_term;
936 for (int t_aux = start; t_aux < end; ++t_aux)
937 {
938 const int t = time_steps - t_aux;
939 const int real_order = std::min(bdf_order, t);
940 const double beta = time_integrator::BDF::betas(real_order - 1);
941
942 Eigen::VectorXd cur_p = adjoint_p.col(t);
943 cur_p(state.boundary_nodes).setZero();
944
945 state.solve_data.damping_form->force_material_derivative(t * dt + t0, state.diff_cached.u(t), state.diff_cached.u(t - 1), -cur_p, damping_term);
946 local_storage.vec += (beta * dt) * damping_term;
947 }
948 });
949
950 for (const LocalThreadVecStorage &local_storage : storage)
951 one_form += local_storage.vec;
952 }
953
954 void AdjointTools::dJ_initial_condition_adjoint_term(
955 const State &state,
956 const Eigen::MatrixXd &adjoint_nu,
957 const Eigen::MatrixXd &adjoint_p,
958 Eigen::VectorXd &one_form)
959 {
960 const int ndof = state.ndof();
961 one_form.setZero(ndof * 2); // half for initial solution, half for initial velocity
962
963 // \partial_q \hat{J}^0 - p_0^T \partial_q g^v - \mu_0^T \partial_q g^u
964 one_form.segment(0, ndof) = -adjoint_nu.col(0); // adjoint_nu[0] actually stores adjoint_mu[0]
965 one_form.segment(ndof, ndof) = -adjoint_p.col(0);
966
967 for (int b : state.boundary_nodes)
968 {
969 one_form(b) = 0;
970 one_form(ndof + b) = 0;
971 }
972 }
973
974 void AdjointTools::dJ_dirichlet_static_adjoint_term(
975 const State &state,
976 const Eigen::MatrixXd &adjoint,
977 Eigen::VectorXd &one_form)
978 {
979 const int n_dirichlet_dof = state.boundary_nodes.size();
980 StiffnessMatrix gradd_h = state.diff_cached.gradu_h(0);
981 std::set<int> boundary_nodes_set(state.boundary_nodes.begin(), state.boundary_nodes.end());
982 gradd_h.prune([&boundary_nodes_set](const Eigen::Index &row, const Eigen::Index &col, const FullNLProblem::Scalar &value) {
983 if (row != col)
984 return value;
985 if (boundary_nodes_set.find(row) == boundary_nodes_set.end())
986 return value;
987 return 0.0;
988 });
989 one_form.setZero(state.ndof());
990 one_form(state.boundary_nodes) -= adjoint(state.boundary_nodes, 0);
991 one_form(state.boundary_nodes) -= (gradd_h.transpose() * adjoint.col(0))(state.boundary_nodes);
992 one_form = utils::flatten(utils::unflatten(one_form, state.mesh->dimension())(state.primitive_to_node(), Eigen::all));
993 }
994
995 void AdjointTools::dJ_dirichlet_transient_adjoint_term(
996 const State &state,
997 const Eigen::MatrixXd &adjoint_nu,
998 const Eigen::MatrixXd &adjoint_p,
999 Eigen::VectorXd &one_form)
1000 {
1001 const double dt = state.args["time"]["dt"];
1002 const int time_steps = state.args["time"]["time_steps"];
1003 const int bdf_order = get_bdf_order(state);
1004 const int n_dirichlet_dof = state.boundary_nodes.size();
1005
1006 // Map dirichlet gradient on each node to dirichlet gradient on boundary ids
1007
1008 one_form.setZero(time_steps * n_dirichlet_dof);
1009 for (int i = 1; i <= time_steps; ++i)
1010 {
1011 const int real_order = std::min(bdf_order, i);
1012 const double beta_dt = time_integrator::BDF::betas(real_order - 1) * dt;
1013
1014 one_form.segment((i - 1) * n_dirichlet_dof, n_dirichlet_dof) = -(1. / beta_dt) * adjoint_p(state.boundary_nodes, i);
1015 }
1016 }
1017
1018 void AdjointTools::dJ_pressure_static_adjoint_term(
1019 const State &state,
1020 const std::vector<int> &boundary_ids,
1021 const Eigen::MatrixXd &sol,
1022 const Eigen::MatrixXd &adjoint,
1023 Eigen::VectorXd &one_form)
1024 {
1025 const int n_pressure_dof = boundary_ids.size();
1026
1027 one_form.setZero(n_pressure_dof);
1028
1029 for (int i = 0; i < boundary_ids.size(); ++i)
1030 {
1031 double pressure_term = state.solve_data.pressure_form->force_pressure_derivative(
1032 state.n_geom_bases,
1033 0,
1034 boundary_ids[i],
1035 sol,
1036 adjoint);
1037 one_form(i) = pressure_term;
1038 }
1039 }
1040
1041 void AdjointTools::dJ_pressure_transient_adjoint_term(
1042 const State &state,
1043 const std::vector<int> &boundary_ids,
1044 const Eigen::MatrixXd &adjoint_nu,
1045 const Eigen::MatrixXd &adjoint_p,
1046 Eigen::VectorXd &one_form)
1047 {
1048 const double t0 = state.args["time"]["t0"];
1049 const double dt = state.args["time"]["dt"];
1050 const int time_steps = state.args["time"]["time_steps"];
1051 const int bdf_order = get_bdf_order(state);
1052
1053 const int n_pressure_dof = boundary_ids.size();
1054
1055 one_form.setZero(time_steps * n_pressure_dof);
1056 Eigen::VectorXd cur_p, cur_nu;
1057 for (int i = time_steps; i > 0; --i)
1058 {
1059 const int real_order = std::min(bdf_order, i);
1060 double beta = time_integrator::BDF::betas(real_order - 1);
1061 double beta_dt = beta * dt;
1062 const double t = i * dt + t0;
1063
1064 cur_p = adjoint_p.col(i);
1065 cur_nu = adjoint_nu.col(i);
1066 cur_p(state.boundary_nodes).setZero();
1067 cur_nu(state.boundary_nodes).setZero();
1068
1069 for (int b = 0; b < boundary_ids.size(); ++b)
1070 {
1071 double pressure_term = state.solve_data.pressure_form->force_pressure_derivative(
1072 state.n_geom_bases,
1073 t,
1074 boundary_ids[b],
1075 state.diff_cached.u(i),
1076 cur_p);
1077 one_form((i - 1) * n_pressure_dof + b) = -beta_dt * pressure_term;
1078 }
1079 }
1080 }
1081
1082 void AdjointTools::dJ_du_step(
1083 const State &state,
1084 const IntegrableFunctional &j,
1085 const Eigen::MatrixXd &solution,
1086 const std::set<int> &interested_ids,
1087 const SpatialIntegralType spatial_integral_type,
1088 const int cur_step,
1089 Eigen::VectorXd &term)
1090 {
1091 const auto &bases = state.bases;
1092 const auto &gbases = state.geom_bases();
1093
1094 const int dim = state.mesh->dimension();
1095 const int actual_dim = state.problem->is_scalar() ? 1 : dim;
1096 const int n_elements = int(bases.size());
1097 const double t0 = state.problem->is_time_dependent() ? state.args["time"]["t0"].get<double>() : 0.0;
1098 const double dt = state.problem->is_time_dependent() ? state.args["time"]["dt"].get<double>() : 0.0;
1099
1100 term = Eigen::MatrixXd::Zero(state.n_bases * actual_dim, 1);
1101
1102 if (!j.depend_on_u() && !j.depend_on_gradu() && !j.depend_on_gradu_local())
1103 return;
1104
1105 if (spatial_integral_type == SpatialIntegralType::Volume)
1106 {
1107 auto storage = utils::create_thread_storage(LocalThreadVecStorage(term.size()));
1108 utils::maybe_parallel_for(n_elements, [&](int start, int end, int thread_id) {
1109 LocalThreadVecStorage &local_storage = utils::get_local_thread_storage(storage, thread_id);
1110
1111 Eigen::MatrixXd u, grad_u;
1112 Eigen::MatrixXd lambda, mu;
1113 Eigen::MatrixXd dj_du, dj_dgradu, dj_dgradx;
1114
1115 IntegrableFunctional::ParameterType params;
1116 params.t = dt * cur_step + t0;
1117 params.step = cur_step;
1118
1119 for (int e = start; e < end; ++e)
1120 {
1121 if (interested_ids.size() != 0 && interested_ids.find(state.mesh->get_body_id(e)) == interested_ids.end())
1122 continue;
1123
1124 assembler::ElementAssemblyValues &vals = local_storage.vals;
1125 state.ass_vals_cache.compute(e, state.mesh->is_volume(), bases[e], gbases[e], vals);
1126
1127 const quadrature::Quadrature &quadrature = vals.quadrature;
1128 local_storage.da = vals.det.array() * quadrature.weights.array();
1129
1130 const Eigen::MatrixXd lame_params = extract_lame_params(state.assembler->parameters(), e, params.t, quadrature.points, vals.val);
1131
1132 const int n_loc_bases_ = int(vals.basis_values.size());
1133
1134 io::Evaluator::interpolate_at_local_vals(e, dim, actual_dim, vals, solution, u, grad_u);
1135
1136 params.elem = e;
1137 params.body_id = state.mesh->get_body_id(e);
1138
1139 dj_dgradu.resize(0, 0);
1140 if (j.depend_on_gradu())
1141 {
1142 j.dj_dgradu(lame_params, quadrature.points, vals.val, u, grad_u, Eigen::MatrixXd::Zero(0, 0) /*Not used*/, vals, params, dj_dgradu);
1143 for (int q = 0; q < dj_dgradu.rows(); q++)
1144 dj_dgradu.row(q) *= local_storage.da(q);
1145 }
1146
1147 dj_du.resize(0, 0);
1148 if (j.depend_on_u())
1149 {
1150 j.dj_du(lame_params, quadrature.points, vals.val, u, grad_u, Eigen::MatrixXd::Zero(0, 0) /*Not used*/, vals, params, dj_du);
1151 for (int q = 0; q < dj_du.rows(); q++)
1152 dj_du.row(q) *= local_storage.da(q);
1153 }
1154
1155 for (int i = 0; i < n_loc_bases_; ++i)
1156 {
1157 const assembler::AssemblyValues &v = vals.basis_values[i];
1158 assert(v.global.size() == 1);
1159 for (int d = 0; d < actual_dim; d++)
1160 {
1161 double val = 0;
1162
1163 // j = j(..., grad u)
1164 if (j.depend_on_gradu())
1165 {
1166 for (int q = 0; q < local_storage.da.size(); ++q)
1167 val += dot(dj_dgradu.block(q, d * dim, 1, dim), v.grad_t_m.row(q));
1168 }
1169
1170 // j = j(..., u)
1171 if (j.depend_on_u())
1172 {
1173 for (int q = 0; q < local_storage.da.size(); ++q)
1174 val += dj_du(q, d) * v.val(q);
1175 }
1176 local_storage.vec(v.global[0].index * actual_dim + d) += val;
1177 }
1178 }
1179 }
1180 });
1181 for (const LocalThreadVecStorage &local_storage : storage)
1182 term += local_storage.vec;
1183 }
1184 else if (spatial_integral_type == SpatialIntegralType::Surface)
1185 {
1186 auto storage = utils::create_thread_storage(LocalThreadVecStorage(term.size()));
1187 utils::maybe_parallel_for(state.total_local_boundary.size(), [&](int start, int end, int thread_id) {
1188 LocalThreadVecStorage &local_storage = utils::get_local_thread_storage(storage, thread_id);
1189
1190 Eigen::MatrixXd uv, samples, gtmp;
1191 Eigen::MatrixXd points, normal;
1192 Eigen::VectorXd weights;
1193
1194 Eigen::MatrixXd u, grad_u;
1195 Eigen::MatrixXd lambda, mu;
1196 Eigen::MatrixXd dj_du, dj_dgradu, dj_dgradu_local;
1197
1198 IntegrableFunctional::ParameterType params;
1199 params.t = dt * cur_step + t0;
1200 params.step = cur_step;
1201
1202 for (int lb_id = start; lb_id < end; ++lb_id)
1203 {
1204 const auto &lb = state.total_local_boundary[lb_id];
1205 const int e = lb.element_id();
1206
1207 for (int i = 0; i < lb.size(); i++)
1208 {
1209 const int global_primitive_id = lb.global_primitive_id(i);
1210 if (interested_ids.size() != 0 && interested_ids.find(state.mesh->get_boundary_id(global_primitive_id)) == interested_ids.end())
1211 continue;
1212
1213 utils::BoundarySampler::boundary_quadrature(lb, state.n_boundary_samples(), *state.mesh, i, false, uv, points, normal, weights);
1214
1215 assembler::ElementAssemblyValues &vals = local_storage.vals;
1216 vals.compute(e, state.mesh->is_volume(), points, bases[e], gbases[e]);
1217 io::Evaluator::interpolate_at_local_vals(e, dim, actual_dim, vals, solution, u, grad_u);
1218
1219 const Eigen::MatrixXd lame_params = extract_lame_params(state.assembler->parameters(), e, params.t, points, vals.val);
1220
1221 // normal = normal * vals.jac_it[0]; // assuming linear geometry
1222
1223 const int n_loc_bases_ = int(vals.basis_values.size());
1224
1225 params.elem = e;
1226 params.body_id = state.mesh->get_body_id(e);
1227 params.boundary_id = state.mesh->get_boundary_id(global_primitive_id);
1228
1229 dj_dgradu.resize(0, 0);
1230 if (j.depend_on_gradu())
1231 {
1232 j.dj_dgradu(lame_params, points, vals.val, u, grad_u, normal, vals, params, dj_dgradu);
1233 for (int q = 0; q < dj_dgradu.rows(); q++)
1234 dj_dgradu.row(q) *= weights(q);
1235 }
1236
1237 dj_dgradu_local.resize(0, 0);
1238 if (j.depend_on_gradu_local())
1239 {
1240 j.dj_dgradu_local(lame_params, points, vals.val, u, grad_u, normal, vals, params, dj_dgradu_local);
1241 for (int q = 0; q < dj_dgradu_local.rows(); q++)
1242 dj_dgradu_local.row(q) *= weights(q);
1243 }
1244
1245 dj_du.resize(0, 0);
1246 if (j.depend_on_u())
1247 {
1248 j.dj_du(lame_params, points, vals.val, u, grad_u, normal, vals, params, dj_du);
1249 for (int q = 0; q < dj_du.rows(); q++)
1250 dj_du.row(q) *= weights(q);
1251 }
1252
1253 for (int l = 0; l < lb.size(); ++l)
1254 {
1255 const auto nodes = bases[e].local_nodes_for_primitive(lb.global_primitive_id(l), *state.mesh);
1256
1257 for (long n = 0; n < nodes.size(); ++n)
1258 {
1259 const assembler::AssemblyValues &v = vals.basis_values[nodes(n)];
1260 assert(v.global.size() == 1);
1261 for (int d = 0; d < actual_dim; d++)
1262 {
1263 double val = 0;
1264
1265 // j = j(x, grad u)
1266 if (j.depend_on_gradu())
1267 {
1268 for (int q = 0; q < weights.size(); ++q)
1269 val += dot(dj_dgradu.block(q, d * dim, 1, dim), v.grad_t_m.row(q));
1270 }
1271 // j = j(x, grad u)
1272 if (j.depend_on_gradu_local())
1273 {
1274 for (int q = 0; q < weights.size(); ++q)
1275 val += dot(dj_dgradu_local.block(q, d * dim, 1, dim), v.grad.row(q));
1276 }
1277 // j = j(x, u)
1278 if (j.depend_on_u())
1279 {
1280 for (int q = 0; q < weights.size(); ++q)
1281 val += dj_du(q, d) * v.val(q);
1282 }
1283 local_storage.vec(v.global[0].index * actual_dim + d) += val;
1284 }
1285 }
1286 }
1287 }
1288 }
1289 });
1290 for (const LocalThreadVecStorage &local_storage : storage)
1291 term += local_storage.vec;
1292 }
1293 else if (spatial_integral_type == SpatialIntegralType::VertexSum)
1294 {
1295 std::vector<bool> traversed(state.n_bases, false);
1296 IntegrableFunctional::ParameterType params;
1297 params.t = dt * cur_step + t0;
1298 params.step = cur_step;
1299 for (int e = 0; e < bases.size(); e++)
1300 {
1301 const auto &bs = bases[e];
1302 for (int i = 0; i < bs.bases.size(); i++)
1303 {
1304 const auto &b = bs.bases[i];
1305 assert(b.global().size() == 1);
1306 const auto &g = b.global()[0];
1307 if (traversed[g.index])
1308 continue;
1309
1310 const Eigen::MatrixXd lame_params = extract_lame_params(state.assembler->parameters(), e, params.t, Eigen::MatrixXd::Zero(1, dim) /*Not used*/, g.node);
1311
1312 params.node = g.index;
1313 params.elem = e;
1314 params.body_id = state.mesh->get_body_id(e);
1315 Eigen::MatrixXd val;
1316 j.dj_du(lame_params, Eigen::MatrixXd::Zero(1, dim) /*Not used*/, g.node, solution.block(g.index * dim, 0, dim, 1).transpose(), Eigen::MatrixXd::Zero(1, dim * actual_dim) /*Not used*/, Eigen::MatrixXd::Zero(0, 0) /*Not used*/, assembler::ElementAssemblyValues(), params, val);
1317 term.block(g.index * actual_dim, 0, actual_dim, 1) += val.transpose();
1318 traversed[g.index] = true;
1319 }
1320 }
1321 }
1322 }
1323
1324 Eigen::VectorXd AdjointTools::map_primitive_to_node_order(const State &state, const Eigen::VectorXd &primitives)
1325 {
1326 int dim = state.mesh->dimension();
1327 assert(primitives.size() == (state.n_geom_bases * dim));
1328 Eigen::VectorXd nodes(primitives.size());
1329 auto map = state.primitive_to_node();
1330 for (int v = 0; v < state.n_geom_bases; ++v)
1331 nodes.segment(map[v] * dim, dim) = primitives.segment(v * dim, dim);
1332 return nodes;
1333 }
1334
1335 Eigen::VectorXd AdjointTools::map_node_to_primitive_order(const State &state, const Eigen::VectorXd &nodes)
1336 {
1337 int dim = state.mesh->dimension();
1338 assert(nodes.size() == (state.n_geom_bases * dim));
1339 Eigen::VectorXd primitives(nodes.size());
1340 auto map = state.node_to_primitive();
1341 for (int v = 0; v < state.n_geom_bases; ++v)
1342 primitives.segment(map[v] * dim, dim) = nodes.segment(v * dim, dim);
1343 return primitives;
1344 }
1345
1346 Eigen::MatrixXd AdjointTools::edge_normal_gradient(const Eigen::MatrixXd &V)
1347 {
1348 DiffScalarBase::setVariableCount(4);
1349 Eigen::Matrix<Diff, 4, 1> full_diff(4, 1);
1350 for (int i = 0; i < 2; i++)
1351 for (int j = 0; j < 2; j++)
1352 full_diff(i * 2 + j) = Diff(i * 2 + j, V(i, j));
1353 auto reduced_diff = edge_normal(full_diff);
1354
1355 Eigen::MatrixXd grad(2, 4);
1356 for (int i = 0; i < 2; ++i)
1357 grad.row(i) = reduced_diff[i].getGradient();
1358
1359 return grad;
1360 }
1361
1362 Eigen::MatrixXd AdjointTools::face_normal_gradient(const Eigen::MatrixXd &V)
1363 {
1364 DiffScalarBase::setVariableCount(9);
1365 Eigen::Matrix<Diff, 9, 1> full_diff(9, 1);
1366 for (int i = 0; i < 3; i++)
1367 for (int j = 0; j < 3; j++)
1368 full_diff(i * 3 + j) = Diff(i * 3 + j, V(i, j));
1369 auto reduced_diff = face_normal(full_diff);
1370
1371 Eigen::MatrixXd grad(3, 9);
1372 for (int i = 0; i < 3; ++i)
1373 grad.row(i) = reduced_diff[i].getGradient();
1374
1375 return grad;
1376 }
1377
1378 Eigen::MatrixXd AdjointTools::edge_velocity_divergence(const Eigen::MatrixXd &V)
1379 {
1380 return line_length_grad(V) / line_length<double>(V);
1381 }
1382
1383 Eigen::MatrixXd AdjointTools::face_velocity_divergence(const Eigen::MatrixXd &V)
1384 {
1385 return triangle_area_grad(V) / triangle_area<double>(V);
1386 }
1387
1388 double AdjointTools::triangle_jacobian(const Eigen::VectorXd &v1, const Eigen::VectorXd &v2, const Eigen::VectorXd &v3)
1389 {
1390 Eigen::VectorXd a = v2 - v1, b = v3 - v1;
1391 return a(0) * b(1) - b(0) * a(1);
1392 }
1393
1394 double AdjointTools::tet_determinant(const Eigen::VectorXd &v1, const Eigen::VectorXd &v2, const Eigen::VectorXd &v3, const Eigen::VectorXd &v4)
1395 {
1396 Eigen::Matrix3d mat;
1397 mat.col(0) << v2 - v1;
1398 mat.col(1) << v3 - v1;
1399 mat.col(2) << v4 - v1;
1400 return mat.determinant();
1401 }
1402
1403 void AdjointTools::scaled_jacobian(const Eigen::MatrixXd &V, const Eigen::MatrixXi &F, Eigen::VectorXd &quality)
1404 {
1405 const int dim = F.cols() - 1;
1406
1407 quality.setZero(F.rows());
1408 if (dim == 2)
1409 {
1410 for (int i = 0; i < F.rows(); i++)
1411 {
1412 Eigen::RowVector3d e0;
1413 e0(2) = 0;
1414 e0.head(2) = V.row(F(i, 2)) - V.row(F(i, 1));
1415 Eigen::RowVector3d e1;
1416 e1(2) = 0;
1417 e1.head(2) = V.row(F(i, 0)) - V.row(F(i, 2));
1418 Eigen::RowVector3d e2;
1419 e2(2) = 0;
1420 e2.head(2) = V.row(F(i, 1)) - V.row(F(i, 0));
1421
1422 double l0 = e0.norm();
1423 double l1 = e1.norm();
1424 double l2 = e2.norm();
1425
1426 double A = 0.5 * (e0.cross(e1)).norm();
1427 double Lmax = std::max(l0 * l1, std::max(l1 * l2, l0 * l2));
1428
1429 quality(i) = 2 * A * (2 / sqrt(3)) / Lmax;
1430 }
1431 }
1432 else
1433 {
1434 for (int i = 0; i < F.rows(); i++)
1435 {
1436 Eigen::RowVector3d e0 = V.row(F(i, 1)) - V.row(F(i, 0));
1437 Eigen::RowVector3d e1 = V.row(F(i, 2)) - V.row(F(i, 1));
1438 Eigen::RowVector3d e2 = V.row(F(i, 0)) - V.row(F(i, 2));
1439 Eigen::RowVector3d e3 = V.row(F(i, 3)) - V.row(F(i, 0));
1440 Eigen::RowVector3d e4 = V.row(F(i, 3)) - V.row(F(i, 1));
1441 Eigen::RowVector3d e5 = V.row(F(i, 3)) - V.row(F(i, 2));
1442
1443 double l0 = e0.norm();
1444 double l1 = e1.norm();
1445 double l2 = e2.norm();
1446 double l3 = e3.norm();
1447 double l4 = e4.norm();
1448 double l5 = e5.norm();
1449
1450 double J = std::abs((e0.cross(e3)).dot(e2));
1451
1452 double a1 = l0 * l2 * l3;
1453 double a2 = l0 * l1 * l4;
1454 double a3 = l1 * l2 * l5;
1455 double a4 = l3 * l4 * l5;
1456
1457 double a = std::max({a1, a2, a3, a4, J});
1458 quality(i) = J * sqrt(2) / a;
1459 }
1460 }
1461 }
1462
1463 bool AdjointTools::is_flipped(const Eigen::MatrixXd &V, const Eigen::MatrixXi &F)
1464 {
1465 if (F.cols() == 3)
1466 {
1467 for (int i = 0; i < F.rows(); i++)
1468 if (triangle_jacobian(V.row(F(i, 0)), V.row(F(i, 1)), V.row(F(i, 2))) <= 0)
1469 return true;
1470 }
1471 else if (F.cols() == 4)
1472 {
1473 for (int i = 0; i < F.rows(); i++)
1474 if (tet_determinant(V.row(F(i, 0)), V.row(F(i, 1)), V.row(F(i, 2)), V.row(F(i, 3))) <= 0)
1475 return true;
1476 }
1477 else
1478 {
1479 return true;
1480 }
1481
1482 return false;
1483 }
1484} // namespace polyfem::solver
V
int V
Definition AdjointNLProblem.cpp:46
vec
Eigen::MatrixXd vec
Definition AdjointTools.cpp:76
val
double val
Definition AdjointTools.cpp:63
vals
assembler::ElementAssemblyValues vals
Definition AdjointTools.cpp:64
vec
Eigen::MatrixXd vec
Definition Assembler.cpp:72
val
double val
Definition Assembler.cpp:86
da
QuadratureVector da
Definition Assembler.cpp:23
vals
ElementAssemblyValues vals
Definition Assembler.cpp:22
gvals
assembler::ElementAssemblyValues gvals
Definition BodyForm.cpp:25
BodyForm.hpp
quadrature
Quadrature quadrature
Definition MassMatrixAssembler.cpp:30
polyfem::IntegrableFunctional::dj_du
void dj_du(const Eigen::MatrixXd &elastic_params, const Eigen::MatrixXd &local_pts, const Eigen::MatrixXd &pts, const Eigen::MatrixXd &u, const Eigen::MatrixXd &grad_u, const Eigen::MatrixXd &reference_normals, const assembler::ElementAssemblyValues &vals, ParameterType &params, Eigen::MatrixXd &val) const
Definition IntegrableFunctional.cpp:50
polyfem::IntegrableFunctional::dj_dgradu
void dj_dgradu(const Eigen::MatrixXd &elastic_params, const Eigen::MatrixXd &local_pts, const Eigen::MatrixXd &pts, const Eigen::MatrixXd &u, const Eigen::MatrixXd &grad_u, const Eigen::MatrixXd &reference_normals, const assembler::ElementAssemblyValues &vals, ParameterType &params, Eigen::MatrixXd &val) const
Definition IntegrableFunctional.cpp:56
polyfem::IntegrableFunctional::evaluate
void evaluate(const Eigen::MatrixXd &elastic_params, const Eigen::MatrixXd &local_pts, const Eigen::MatrixXd &pts, const Eigen::MatrixXd &u, const Eigen::MatrixXd &grad_u, const Eigen::MatrixXd &reference_normals, const assembler::ElementAssemblyValues &vals, ParameterType &params, Eigen::MatrixXd &val) const
Definition IntegrableFunctional.cpp:38
polyfem::IntegrableFunctional::dj_dx
void dj_dx(const Eigen::MatrixXd &elastic_params, const Eigen::MatrixXd &local_pts, const Eigen::MatrixXd &pts, const Eigen::MatrixXd &u, const Eigen::MatrixXd &grad_u, const Eigen::MatrixXd &reference_normals, const assembler::ElementAssemblyValues &vals, ParameterType &params, Eigen::MatrixXd &val) const
Definition IntegrableFunctional.cpp:44
polyfem::State
main class that contains the polyfem solver and all its state
Definition State.hpp:79
polyfem::State::initial_vel_update
Eigen::MatrixXd initial_vel_update
Definition State.hpp:704
polyfem::State::n_bases
int n_bases
number of bases
Definition State.hpp:178
polyfem::State::ass_vals_cache
assembler::AssemblyValsCache ass_vals_cache
used to store assembly values for small problems
Definition State.hpp:196
polyfem::State::is_adhesion_enabled
bool is_adhesion_enabled() const
does the simulation have adhesion
Definition State.hpp:564
polyfem::State::geom_bases
const std::vector< basis::ElementBases > & geom_bases() const
Get a constant reference to the geometry mapping bases.
Definition State.hpp:223
polyfem::State::assembler
std::shared_ptr< assembler::Assembler > assembler
assemblers
Definition State.hpp:155
polyfem::State::collision_mesh
ipc::CollisionMesh collision_mesh
IPC collision mesh.
Definition State.hpp:520
polyfem::State::mass_matrix_assembler
std::shared_ptr< assembler::Mass > mass_matrix_assembler
Definition State.hpp:157
polyfem::State::primitive_to_node
std::vector< int > primitive_to_node() const
Definition State.cpp:227
polyfem::State::mesh
std::unique_ptr< mesh::Mesh > mesh
current mesh, it can be a Mesh2D or Mesh3D
Definition State.hpp:471
polyfem::State::mesh_nodes
std::shared_ptr< polyfem::mesh::MeshNodes > mesh_nodes
Mapping from input nodes to FE nodes.
Definition State.hpp:193
polyfem::State::problem
std::shared_ptr< assembler::Problem > problem
current problem, it contains rhs and bc
Definition State.hpp:168
polyfem::State::node_to_primitive
std::vector< int > node_to_primitive() const
Definition State.cpp:234
polyfem::State::args
json args
main input arguments containing all defaults
Definition State.hpp:101
polyfem::State::diff_cached
solver::DiffCache diff_cached
Definition State.hpp:656
polyfem::State::initial_velocity
void initial_velocity(Eigen::MatrixXd &velocity) const
Load or compute the initial velocity.
Definition StateSolve.cpp:90
polyfem::State::initial_acceleration
void initial_acceleration(Eigen::MatrixXd &acceleration) const
Load or compute the initial acceleration.
Definition StateSolve.cpp:103
polyfem::State::bases
std::vector< basis::ElementBases > bases
FE bases, the size is #elements.
Definition State.hpp:171
polyfem::State::n_geom_bases
int n_geom_bases
number of geometric bases
Definition State.hpp:182
polyfem::State::total_local_boundary
std::vector< mesh::LocalBoundary > total_local_boundary
mapping from elements to nodes for all mesh
Definition State.hpp:436
polyfem::State::ndof
int ndof() const
Definition State.hpp:660
polyfem::State::mass_ass_vals_cache
assembler::AssemblyValsCache mass_ass_vals_cache
Definition State.hpp:197
polyfem::State::boundary_nodes
std::vector< int > boundary_nodes
list of boundary nodes
Definition State.hpp:432
polyfem::State::solve_data
solver::SolveData solve_data
timedependent stuff cached
Definition State.hpp:321
polyfem::State::is_contact_enabled
bool is_contact_enabled() const
does the simulation have contact
Definition State.hpp:556
polyfem::State::basis_nodes_to_gbasis_nodes
StiffnessMatrix basis_nodes_to_gbasis_nodes
Definition State.hpp:706
polyfem::assembler::AssemblyValsCache::compute
void compute(const int el_index, const bool is_volume, const basis::ElementBases &basis, const basis::ElementBases &gbasis, ElementAssemblyValues &vals) const
retrieves cached basis evaluation and geometric for the given element if it doesn't exist,...
Definition AssemblyValsCache.cpp:45
polyfem::assembler::AssemblyValues
stores per local bases evaluations
Definition AssemblyValues.hpp:14
polyfem::assembler::AssemblyValues::global
std::vector< basis::Local2Global > global
Definition AssemblyValues.hpp:20
polyfem::assembler::ElementAssemblyValues
stores per element basis values at given quadrature points and geometric mapping
Definition ElementAssemblyValues.hpp:14
polyfem::assembler::ElementAssemblyValues::compute
void compute(const int el_index, const bool is_volume, const Eigen::MatrixXd &pts, const basis::ElementBases &basis, const basis::ElementBases &gbasis)
computes the per element values at the local (ref el) points (pts) sets basis_values,...
Definition ElementAssemblyValues.cpp:164
polyfem::io::Evaluator::interpolate_at_local_vals
static void interpolate_at_local_vals(const mesh::Mesh &mesh, const bool is_problem_scalar, const std::vector< basis::ElementBases > &bases, const std::vector< basis::ElementBases > &gbases, const int el_index, const Eigen::MatrixXd &local_pts, const Eigen::MatrixXd &fun, Eigen::MatrixXd &result, Eigen::MatrixXd &result_grad)
interpolate solution and gradient at element (calls interpolate_at_local_vals with sol)
Definition Evaluator.cpp:571
polyfem::solver::DiffCache::collision_set
const ipc::NormalCollisions & collision_set(int step) const
Definition DiffCache.hpp:158
polyfem::solver::DiffCache::gradu_h
const StiffnessMatrix & gradu_h(int step) const
Definition DiffCache.hpp:149
polyfem::solver::DiffCache::friction_collision_set
const ipc::TangentialCollisions & friction_collision_set(int step) const
Definition DiffCache.hpp:172
polyfem::solver::DiffCache::disp_grad
Eigen::MatrixXd disp_grad(int step=0) const
Definition DiffCache.hpp:125
polyfem::solver::DiffCache::smooth_collision_set
const ipc::SmoothCollisions & smooth_collision_set(int step) const
Definition DiffCache.hpp:165
polyfem::solver::DiffCache::tangential_adhesion_collision_set
const ipc::TangentialCollisions & tangential_adhesion_collision_set(int step) const
Definition DiffCache.hpp:186
polyfem::solver::DiffCache::normal_adhesion_collision_set
const ipc::NormalCollisions & normal_adhesion_collision_set(int step) const
Definition DiffCache.hpp:179
polyfem::solver::DiffCache::v
Eigen::VectorXd v(int step) const
Definition DiffCache.hpp:134
polyfem::solver::DiffCache::u
Eigen::VectorXd u(int step) const
Definition DiffCache.hpp:127
polyfem::solver::PeriodicMeshToMesh::apply_jacobian
Eigen::VectorXd apply_jacobian(const Eigen::VectorXd &grad, const Eigen::VectorXd &x) const override
Definition PeriodicMeshToMesh.cpp:137
polyfem::solver::SolveData::friction_form
std::shared_ptr< solver::FrictionForm > friction_form
Definition SolveData.hpp:178
polyfem::solver::SolveData::inertia_form
std::shared_ptr< solver::InertiaForm > inertia_form
Definition SolveData.hpp:179
polyfem::solver::SolveData::periodic_contact_form
std::shared_ptr< solver::PeriodicContactForm > periodic_contact_form
Definition SolveData.hpp:184
polyfem::solver::SolveData::pressure_form
std::shared_ptr< solver::PressureForm > pressure_form
Definition SolveData.hpp:180
polyfem::solver::SolveData::body_form
std::shared_ptr< solver::BodyForm > body_form
Definition SolveData.hpp:174
polyfem::solver::SolveData::nl_problem
std::shared_ptr< solver::NLProblem > nl_problem
Definition SolveData.hpp:170
polyfem::solver::SolveData::normal_adhesion_form
std::shared_ptr< solver::NormalAdhesionForm > normal_adhesion_form
Definition SolveData.hpp:181
polyfem::solver::SolveData::contact_form
std::shared_ptr< solver::ContactForm > contact_form
Definition SolveData.hpp:175
polyfem::solver::SolveData::damping_form
std::shared_ptr< solver::ElasticForm > damping_form
Definition SolveData.hpp:176
polyfem::solver::SolveData::elastic_form
std::shared_ptr< solver::ElasticForm > elastic_form
Definition SolveData.hpp:177
polyfem::solver::SolveData::tangential_adhesion_form
std::shared_ptr< solver::TangentialAdhesionForm > tangential_adhesion_form
Definition SolveData.hpp:182
generate_mooney_rivlin.grad
grad
Definition generate_mooney_rivlin.py:51
generate_rotation_mtx.T
T
Definition generate_rotation_mtx.py:43
polyfem::assembler::cross
Eigen::Matrix< double, dim, 1 > cross(const Eigen::Matrix< double, dim, 1 > &x, const Eigen::Matrix< double, dim, 1 > &y)
Definition NeoHookeanElasticity.cpp:439
polyfem::solver::AdjointTools::tet_determinant
double tet_determinant(const Eigen::VectorXd &v1, const Eigen::VectorXd &v2, const Eigen::VectorXd &v3, const Eigen::VectorXd &v4)
Definition AdjointTools.cpp:1394
polyfem::solver::AdjointTools::integrate_objective
double integrate_objective(const State &state, const IntegrableFunctional &j, const Eigen::MatrixXd &solution, const std::set< int > &interested_ids, const SpatialIntegralType spatial_integral_type, const int cur_step=0)
Definition AdjointTools.cpp:182
polyfem::solver::AdjointTools::dJ_shape_homogenization_adjoint_term
void dJ_shape_homogenization_adjoint_term(const State &state, const Eigen::MatrixXd &sol, const Eigen::MatrixXd &adjoint, Eigen::VectorXd &one_form)
Definition AdjointTools.cpp:610
polyfem::solver::AdjointTools::triangle_jacobian
double triangle_jacobian(const Eigen::VectorXd &v1, const Eigen::VectorXd &v2, const Eigen::VectorXd &v3)
Definition AdjointTools.cpp:1388
polyfem::solver::Diff
DScalar1< double, Eigen::Matrix< double, Eigen::Dynamic, 1 > > Diff
Definition SurfaceTractionForms.cpp:236
polyfem::utils::vector2matrix
void vector2matrix(const Eigen::VectorXd &vec, Eigen::MatrixXd &mat)
Definition MatrixUtils.cpp:69
polyfem::utils::get_local_thread_storage
auto & get_local_thread_storage(Storages &storage, int thread_id)
polyfem::utils::create_thread_storage
auto create_thread_storage(const LocalStorage &initial_local_storage)
polyfem::utils::triangle_area
double triangle_area(const Eigen::MatrixXd V)
Compute the signed area of a triangle defined by three points.
Definition GeometryUtils.cpp:30
polyfem::utils::maybe_parallel_for
void maybe_parallel_for(int size, const std::function< void(int, int, int)> &partial_for)
polyfem::QuadratureVector
Eigen::Matrix< double, Eigen::Dynamic, 1, 0, MAX_QUAD_POINTS, 1 > QuadratureVector
Definition Types.hpp:17
polyfem::log_and_throw_adjoint_error
void log_and_throw_adjoint_error(const std::string &msg)
Definition Logger.cpp:79
polyfem::StiffnessMatrix
Eigen::SparseMatrix< double, Eigen::ColMajor > StiffnessMatrix
Definition Types.hpp:22
DScalar1
Automatic differentiation scalar with first-order derivatives.
Definition autodiff.h:112
DiffScalarBase::setVariableCount
static void setVariableCount(size_t value)
Set the independent variable count used by the automatic differentiation layer.
Definition autodiff.h:54
polyfem::IntegrableFunctional::ParameterType
Parameters for the functional evaluation.
Definition IntegrableFunctional.hpp:18