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