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 ipc::BarrierPotential bp = barrier_contact->barrier_potential();
936 bp.set_stiffness(barrier_contact->barrier_stiffness());
937 Eigen::MatrixXd force = state.collision_mesh.to_full_dof(
938 -state.solve_data.friction_form->friction_potential().force(
939 diff_cache.friction_collision_set(t),
940 state.collision_mesh,
941 state.collision_mesh.rest_positions(),
942 /*lagged_displacements=*/surface_solution_prev,
943 surface_velocities,
944 bp,
945 0., true));
946
947 Eigen::VectorXd cur_p = adjoint_p.col(t);
948 cur_p(state.boundary_nodes).setZero();
949
950 one_form(0) += dot(cur_p, force) * beta * dt;
951 }
952 }
953 }
954
955 void AdjointTools::dJ_damping_transient_adjoint_term(
956 const State &state,
957 const DiffCache &diff_cache,
958 const Eigen::MatrixXd &adjoint_nu,
959 const Eigen::MatrixXd &adjoint_p,
960 Eigen::VectorXd &one_form)
961 {
962 const double t0 = state.args["time"]["t0"];
963 const double dt = state.args["time"]["dt"];
964 const int time_steps = state.args["time"]["time_steps"];
965 const int bdf_order = get_bdf_order(state);
966
967 one_form.setZero(2);
968
969 auto storage = utils::create_thread_storage(LocalThreadVecStorage(one_form.size()));
970
971 utils::maybe_parallel_for(time_steps, [&](int start, int end, int thread_id) {
972 LocalThreadVecStorage &local_storage = utils::get_local_thread_storage(storage, thread_id);
973 Eigen::VectorXd damping_term;
974 for (int t_aux = start; t_aux < end; ++t_aux)
975 {
976 const int t = time_steps - t_aux;
977 const int real_order = std::min(bdf_order, t);
978 const double beta = time_integrator::BDF::betas(real_order - 1);
979
980 Eigen::VectorXd cur_p = adjoint_p.col(t);
981 cur_p(state.boundary_nodes).setZero();
982
983 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);
984 local_storage.vec += (beta * dt) * damping_term;
985 }
986 });
987
988 for (const LocalThreadVecStorage &local_storage : storage)
989 one_form += local_storage.vec;
990 }
991
992 void AdjointTools::dJ_initial_condition_adjoint_term(
993 const State &state,
994 const Eigen::MatrixXd &adjoint_nu,
995 const Eigen::MatrixXd &adjoint_p,
996 Eigen::VectorXd &one_form)
997 {
998 const int ndof = state.ndof();
999 one_form.setZero(ndof * 2); // half for initial solution, half for initial velocity
1000
1001 // \partial_q \hat{J}^0 - p_0^T \partial_q g^v - \mu_0^T \partial_q g^u
1002 one_form.segment(0, ndof) = -adjoint_nu.col(0); // adjoint_nu[0] actually stores adjoint_mu[0]
1003 one_form.segment(ndof, ndof) = -adjoint_p.col(0);
1004
1005 for (int b : state.boundary_nodes)
1006 {
1007 one_form(b) = 0;
1008 one_form(ndof + b) = 0;
1009 }
1010 }
1011
1012 void AdjointTools::dJ_dirichlet_static_adjoint_term(
1013 const State &state,
1014 const DiffCache &diff_cache,
1015 const Eigen::MatrixXd &adjoint,
1016 Eigen::VectorXd &one_form)
1017 {
1018 const int n_dirichlet_dof = state.boundary_nodes.size();
1019 StiffnessMatrix gradd_h = diff_cache.gradu_h(0);
1020 std::set<int> boundary_nodes_set(state.boundary_nodes.begin(), state.boundary_nodes.end());
1021 gradd_h.prune([&boundary_nodes_set](const Eigen::Index &row, const Eigen::Index &col, const FullNLProblem::Scalar &value) {
1022 if (row != col)
1023 return value;
1024 if (boundary_nodes_set.find(row) == boundary_nodes_set.end())
1025 return value;
1026 return 0.0;
1027 });
1028 one_form.setZero(state.ndof());
1029 one_form(state.boundary_nodes) -= adjoint(state.boundary_nodes, 0);
1030 one_form(state.boundary_nodes) -= (gradd_h.transpose() * adjoint.col(0))(state.boundary_nodes);
1031 one_form = utils::flatten(utils::unflatten(one_form, state.mesh->dimension())(state.primitive_to_node(), Eigen::all));
1032 }
1033
1034 void AdjointTools::dJ_dirichlet_transient_adjoint_term(
1035 const State &state,
1036 const Eigen::MatrixXd &adjoint_nu,
1037 const Eigen::MatrixXd &adjoint_p,
1038 Eigen::VectorXd &one_form)
1039 {
1040 const double dt = state.args["time"]["dt"];
1041 const int time_steps = state.args["time"]["time_steps"];
1042 const int bdf_order = get_bdf_order(state);
1043 const int n_dirichlet_dof = state.boundary_nodes.size();
1044
1045 // Map dirichlet gradient on each node to dirichlet gradient on boundary ids
1046
1047 one_form.setZero(time_steps * n_dirichlet_dof);
1048 for (int i = 1; i <= time_steps; ++i)
1049 {
1050 const int real_order = std::min(bdf_order, i);
1051 const double beta_dt = time_integrator::BDF::betas(real_order - 1) * dt;
1052
1053 one_form.segment((i - 1) * n_dirichlet_dof, n_dirichlet_dof) = -(1. / beta_dt) * adjoint_p(state.boundary_nodes, i);
1054 }
1055 }
1056
1057 void AdjointTools::dJ_pressure_static_adjoint_term(
1058 const State &state,
1059 const std::vector<int> &boundary_ids,
1060 const Eigen::MatrixXd &sol,
1061 const Eigen::MatrixXd &adjoint,
1062 Eigen::VectorXd &one_form)
1063 {
1064 const int n_pressure_dof = boundary_ids.size();
1065
1066 one_form.setZero(n_pressure_dof);
1067
1068 for (int i = 0; i < boundary_ids.size(); ++i)
1069 {
1070 double pressure_term = PressureForceDerivative::force_pressure_derivative(
1071 *state.solve_data.pressure_form,
1072 state.n_geom_bases,
1073 0,
1074 boundary_ids[i],
1075 sol,
1076 adjoint);
1077 one_form(i) = pressure_term;
1078 }
1079 }
1080
1081 void AdjointTools::dJ_pressure_transient_adjoint_term(
1082 const State &state,
1083 const DiffCache &diff_cache,
1084 const std::vector<int> &boundary_ids,
1085 const Eigen::MatrixXd &adjoint_nu,
1086 const Eigen::MatrixXd &adjoint_p,
1087 Eigen::VectorXd &one_form)
1088 {
1089 const double t0 = state.args["time"]["t0"];
1090 const double dt = state.args["time"]["dt"];
1091 const int time_steps = state.args["time"]["time_steps"];
1092 const int bdf_order = get_bdf_order(state);
1093
1094 const int n_pressure_dof = boundary_ids.size();
1095
1096 one_form.setZero(time_steps * n_pressure_dof);
1097 Eigen::VectorXd cur_p, cur_nu;
1098 for (int i = time_steps; i > 0; --i)
1099 {
1100 const int real_order = std::min(bdf_order, i);
1101 double beta = time_integrator::BDF::betas(real_order - 1);
1102 double beta_dt = beta * dt;
1103 const double t = i * dt + t0;
1104
1105 cur_p = adjoint_p.col(i);
1106 cur_nu = adjoint_nu.col(i);
1107 cur_p(state.boundary_nodes).setZero();
1108 cur_nu(state.boundary_nodes).setZero();
1109
1110 for (int b = 0; b < boundary_ids.size(); ++b)
1111 {
1112 double pressure_term = PressureForceDerivative::force_pressure_derivative(
1113 *state.solve_data.pressure_form,
1114 state.n_geom_bases,
1115 t,
1116 boundary_ids[b],
1117 diff_cache.u(i),
1118 cur_p);
1119 one_form((i - 1) * n_pressure_dof + b) = -beta_dt * pressure_term;
1120 }
1121 }
1122 }
1123
1124 void AdjointTools::dJ_du_step(
1125 const State &state,
1126 const IntegrableFunctional &j,
1127 const Eigen::MatrixXd &solution,
1128 const std::set<int> &interested_ids,
1129 const SpatialIntegralType spatial_integral_type,
1130 const int cur_step,
1131 Eigen::VectorXd &term)
1132 {
1133 const auto &bases = state.bases;
1134 const auto &gbases = state.geom_bases();
1135
1136 const int dim = state.mesh->dimension();
1137 const int actual_dim = state.problem->is_scalar() ? 1 : dim;
1138 const int n_elements = int(bases.size());
1139 const double t0 = state.problem->is_time_dependent() ? state.args["time"]["t0"].get<double>() : 0.0;
1140 const double dt = state.problem->is_time_dependent() ? state.args["time"]["dt"].get<double>() : 0.0;
1141
1142 term = Eigen::MatrixXd::Zero(state.n_bases * actual_dim, 1);
1143
1144 if (!j.depend_on_u() && !j.depend_on_gradu() && !j.depend_on_gradu_local())
1145 return;
1146
1147 if (spatial_integral_type == SpatialIntegralType::Volume)
1148 {
1149 auto storage = utils::create_thread_storage(LocalThreadVecStorage(term.size()));
1150 utils::maybe_parallel_for(n_elements, [&](int start, int end, int thread_id) {
1151 LocalThreadVecStorage &local_storage = utils::get_local_thread_storage(storage, thread_id);
1152
1153 Eigen::MatrixXd u, grad_u;
1154 Eigen::MatrixXd lambda, mu;
1155 Eigen::MatrixXd dj_du, dj_dgradu, dj_dgradx;
1156
1157 IntegrableFunctional::ParameterType params;
1158 params.t = dt * cur_step + t0;
1159 params.step = cur_step;
1160
1161 for (int e = start; e < end; ++e)
1162 {
1163 if (interested_ids.size() != 0 && interested_ids.find(state.mesh->get_body_id(e)) == interested_ids.end())
1164 continue;
1165
1166 assembler::ElementAssemblyValues &vals = local_storage.vals;
1167 state.ass_vals_cache.compute(e, state.mesh->is_volume(), bases[e], gbases[e], vals);
1168
1169 const quadrature::Quadrature &quadrature = vals.quadrature;
1170 local_storage.da = vals.det.array() * quadrature.weights.array();
1171
1172 const Eigen::MatrixXd lame_params = extract_lame_params(state.assembler->parameters(), e, params.t, quadrature.points, vals.val);
1173
1174 const int n_loc_bases_ = int(vals.basis_values.size());
1175
1176 io::Evaluator::interpolate_at_local_vals(e, dim, actual_dim, vals, solution, u, grad_u);
1177
1178 params.elem = e;
1179 params.body_id = state.mesh->get_body_id(e);
1180
1181 dj_dgradu.resize(0, 0);
1182 if (j.depend_on_gradu())
1183 {
1184 j.dj_dgradu(lame_params, quadrature.points, vals.val, u, grad_u, Eigen::MatrixXd::Zero(0, 0) /*Not used*/, vals, params, dj_dgradu);
1185 for (int q = 0; q < dj_dgradu.rows(); q++)
1186 dj_dgradu.row(q) *= local_storage.da(q);
1187 }
1188
1189 dj_du.resize(0, 0);
1190 if (j.depend_on_u())
1191 {
1192 j.dj_du(lame_params, quadrature.points, vals.val, u, grad_u, Eigen::MatrixXd::Zero(0, 0) /*Not used*/, vals, params, dj_du);
1193 for (int q = 0; q < dj_du.rows(); q++)
1194 dj_du.row(q) *= local_storage.da(q);
1195 }
1196
1197 for (int i = 0; i < n_loc_bases_; ++i)
1198 {
1199 const assembler::AssemblyValues &v = vals.basis_values[i];
1200 assert(v.global.size() == 1);
1201 for (int d = 0; d < actual_dim; d++)
1202 {
1203 double val = 0;
1204
1205 // j = j(..., grad u)
1206 if (j.depend_on_gradu())
1207 {
1208 for (int q = 0; q < local_storage.da.size(); ++q)
1209 val += dot(dj_dgradu.block(q, d * dim, 1, dim), v.grad_t_m.row(q));
1210 }
1211
1212 // j = j(..., u)
1213 if (j.depend_on_u())
1214 {
1215 for (int q = 0; q < local_storage.da.size(); ++q)
1216 val += dj_du(q, d) * v.val(q);
1217 }
1218 local_storage.vec(v.global[0].index * actual_dim + d) += val;
1219 }
1220 }
1221 }
1222 });
1223 for (const LocalThreadVecStorage &local_storage : storage)
1224 term += local_storage.vec;
1225 }
1226 else if (spatial_integral_type == SpatialIntegralType::Surface)
1227 {
1228 auto storage = utils::create_thread_storage(LocalThreadVecStorage(term.size()));
1229 utils::maybe_parallel_for(state.total_local_boundary.size(), [&](int start, int end, int thread_id) {
1230 LocalThreadVecStorage &local_storage = utils::get_local_thread_storage(storage, thread_id);
1231
1232 Eigen::MatrixXd uv, samples, gtmp;
1233 Eigen::MatrixXd points, normal;
1234 Eigen::VectorXd weights;
1235
1236 Eigen::MatrixXd u, grad_u;
1237 Eigen::MatrixXd lambda, mu;
1238 Eigen::MatrixXd dj_du, dj_dgradu, dj_dgradu_local;
1239
1240 IntegrableFunctional::ParameterType params;
1241 params.t = dt * cur_step + t0;
1242 params.step = cur_step;
1243
1244 for (int lb_id = start; lb_id < end; ++lb_id)
1245 {
1246 const auto &lb = state.total_local_boundary[lb_id];
1247 const int e = lb.element_id();
1248
1249 for (int i = 0; i < lb.size(); i++)
1250 {
1251 const int global_primitive_id = lb.global_primitive_id(i);
1252 if (interested_ids.size() != 0 && interested_ids.find(state.mesh->get_boundary_id(global_primitive_id)) == interested_ids.end())
1253 continue;
1254
1255 utils::BoundarySampler::boundary_quadrature(lb, state.n_boundary_samples(), *state.mesh, i, false, uv, points, normal, weights);
1256
1257 assembler::ElementAssemblyValues &vals = local_storage.vals;
1258 vals.compute(e, state.mesh->is_volume(), points, bases[e], gbases[e]);
1259 io::Evaluator::interpolate_at_local_vals(e, dim, actual_dim, vals, solution, u, grad_u);
1260
1261 const Eigen::MatrixXd lame_params = extract_lame_params(state.assembler->parameters(), e, params.t, points, vals.val);
1262
1263 // normal = normal * vals.jac_it[0]; // assuming linear geometry
1264
1265 const int n_loc_bases_ = int(vals.basis_values.size());
1266
1267 params.elem = e;
1268 params.body_id = state.mesh->get_body_id(e);
1269 params.boundary_id = state.mesh->get_boundary_id(global_primitive_id);
1270
1271 dj_dgradu.resize(0, 0);
1272 if (j.depend_on_gradu())
1273 {
1274 j.dj_dgradu(lame_params, points, vals.val, u, grad_u, normal, vals, params, dj_dgradu);
1275 for (int q = 0; q < dj_dgradu.rows(); q++)
1276 dj_dgradu.row(q) *= weights(q);
1277 }
1278
1279 dj_dgradu_local.resize(0, 0);
1280 if (j.depend_on_gradu_local())
1281 {
1282 j.dj_dgradu_local(lame_params, points, vals.val, u, grad_u, normal, vals, params, dj_dgradu_local);
1283 for (int q = 0; q < dj_dgradu_local.rows(); q++)
1284 dj_dgradu_local.row(q) *= weights(q);
1285 }
1286
1287 dj_du.resize(0, 0);
1288 if (j.depend_on_u())
1289 {
1290 j.dj_du(lame_params, points, vals.val, u, grad_u, normal, vals, params, dj_du);
1291 for (int q = 0; q < dj_du.rows(); q++)
1292 dj_du.row(q) *= weights(q);
1293 }
1294
1295 for (int l = 0; l < lb.size(); ++l)
1296 {
1297 const auto nodes = bases[e].local_nodes_for_primitive(lb.global_primitive_id(l), *state.mesh);
1298
1299 for (long n = 0; n < nodes.size(); ++n)
1300 {
1301 const assembler::AssemblyValues &v = vals.basis_values[nodes(n)];
1302 assert(v.global.size() == 1);
1303 for (int d = 0; d < actual_dim; d++)
1304 {
1305 double val = 0;
1306
1307 // j = j(x, grad u)
1308 if (j.depend_on_gradu())
1309 {
1310 for (int q = 0; q < weights.size(); ++q)
1311 val += dot(dj_dgradu.block(q, d * dim, 1, dim), v.grad_t_m.row(q));
1312 }
1313 // j = j(x, grad u)
1314 if (j.depend_on_gradu_local())
1315 {
1316 for (int q = 0; q < weights.size(); ++q)
1317 val += dot(dj_dgradu_local.block(q, d * dim, 1, dim), v.grad.row(q));
1318 }
1319 // j = j(x, u)
1320 if (j.depend_on_u())
1321 {
1322 for (int q = 0; q < weights.size(); ++q)
1323 val += dj_du(q, d) * v.val(q);
1324 }
1325 local_storage.vec(v.global[0].index * actual_dim + d) += val;
1326 }
1327 }
1328 }
1329 }
1330 }
1331 });
1332 for (const LocalThreadVecStorage &local_storage : storage)
1333 term += local_storage.vec;
1334 }
1335 else if (spatial_integral_type == SpatialIntegralType::VertexSum)
1336 {
1337 std::vector<bool> traversed(state.n_bases, false);
1338 IntegrableFunctional::ParameterType params;
1339 params.t = dt * cur_step + t0;
1340 params.step = cur_step;
1341 for (int e = 0; e < bases.size(); e++)
1342 {
1343 const auto &bs = bases[e];
1344 for (int i = 0; i < bs.bases.size(); i++)
1345 {
1346 const auto &b = bs.bases[i];
1347 assert(b.global().size() == 1);
1348 const auto &g = b.global()[0];
1349 if (traversed[g.index])
1350 continue;
1351
1352 const Eigen::MatrixXd lame_params = extract_lame_params(state.assembler->parameters(), e, params.t, Eigen::MatrixXd::Zero(1, dim) /*Not used*/, g.node);
1353
1354 params.node = g.index;
1355 params.elem = e;
1356 params.body_id = state.mesh->get_body_id(e);
1357 Eigen::MatrixXd val;
1358 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);
1359 term.block(g.index * actual_dim, 0, actual_dim, 1) += val.transpose();
1360 traversed[g.index] = true;
1361 }
1362 }
1363 }
1364 }
1365
1366 Eigen::VectorXd AdjointTools::map_primitive_to_node_order(const State &state, const Eigen::VectorXd &primitives)
1367 {
1368 int dim = state.mesh->dimension();
1369 assert(primitives.size() == (state.n_geom_bases * dim));
1370 Eigen::VectorXd nodes(primitives.size());
1371 auto map = state.primitive_to_node();
1372 for (int v = 0; v < state.n_geom_bases; ++v)
1373 nodes.segment(map[v] * dim, dim) = primitives.segment(v * dim, dim);
1374 return nodes;
1375 }
1376
1377 Eigen::VectorXd AdjointTools::map_node_to_primitive_order(const State &state, const Eigen::VectorXd &nodes)
1378 {
1379 int dim = state.mesh->dimension();
1380 assert(nodes.size() == (state.n_geom_bases * dim));
1381 Eigen::VectorXd primitives(nodes.size());
1382 auto map = state.node_to_primitive();
1383 for (int v = 0; v < state.n_geom_bases; ++v)
1384 primitives.segment(map[v] * dim, dim) = nodes.segment(v * dim, dim);
1385 return primitives;
1386 }
1387
1388 Eigen::MatrixXd AdjointTools::edge_normal_gradient(const Eigen::MatrixXd &V)
1389 {
1390 DiffScalarBase::setVariableCount(4);
1391 Eigen::Matrix<Diff, 4, 1> full_diff(4, 1);
1392 for (int i = 0; i < 2; i++)
1393 for (int j = 0; j < 2; j++)
1394 full_diff(i * 2 + j) = Diff(i * 2 + j, V(i, j));
1395 auto reduced_diff = edge_normal(full_diff);
1396
1397 Eigen::MatrixXd grad(2, 4);
1398 for (int i = 0; i < 2; ++i)
1399 grad.row(i) = reduced_diff[i].getGradient();
1400
1401 return grad;
1402 }
1403
1404 Eigen::MatrixXd AdjointTools::face_normal_gradient(const Eigen::MatrixXd &V)
1405 {
1406 DiffScalarBase::setVariableCount(9);
1407 Eigen::Matrix<Diff, 9, 1> full_diff(9, 1);
1408 for (int i = 0; i < 3; i++)
1409 for (int j = 0; j < 3; j++)
1410 full_diff(i * 3 + j) = Diff(i * 3 + j, V(i, j));
1411 auto reduced_diff = face_normal(full_diff);
1412
1413 Eigen::MatrixXd grad(3, 9);
1414 for (int i = 0; i < 3; ++i)
1415 grad.row(i) = reduced_diff[i].getGradient();
1416
1417 return grad;
1418 }
1419
1420 Eigen::MatrixXd AdjointTools::edge_velocity_divergence(const Eigen::MatrixXd &V)
1421 {
1422 return line_length_grad(V) / line_length<double>(V);
1423 }
1424
1425 Eigen::MatrixXd AdjointTools::face_velocity_divergence(const Eigen::MatrixXd &V)
1426 {
1427 return triangle_area_grad(V) / triangle_area<double>(V);
1428 }
1429
1430 void AdjointTools::scaled_jacobian(const Eigen::MatrixXd &V, const Eigen::MatrixXi &F, Eigen::VectorXd &quality)
1431 {
1432 const int dim = F.cols() - 1;
1433
1434 quality.setZero(F.rows());
1435 if (dim == 2)
1436 {
1437 for (int i = 0; i < F.rows(); i++)
1438 {
1439 Eigen::RowVector3d e0;
1440 e0(2) = 0;
1441 e0.head(2) = V.row(F(i, 2)) - V.row(F(i, 1));
1442 Eigen::RowVector3d e1;
1443 e1(2) = 0;
1444 e1.head(2) = V.row(F(i, 0)) - V.row(F(i, 2));
1445 Eigen::RowVector3d e2;
1446 e2(2) = 0;
1447 e2.head(2) = V.row(F(i, 1)) - V.row(F(i, 0));
1448
1449 double l0 = e0.norm();
1450 double l1 = e1.norm();
1451 double l2 = e2.norm();
1452
1453 double A = 0.5 * (e0.cross(e1)).norm();
1454 double Lmax = std::max(l0 * l1, std::max(l1 * l2, l0 * l2));
1455
1456 quality(i) = 2 * A * (2 / sqrt(3)) / Lmax;
1457 }
1458 }
1459 else
1460 {
1461 for (int i = 0; i < F.rows(); i++)
1462 {
1463 Eigen::RowVector3d e0 = V.row(F(i, 1)) - V.row(F(i, 0));
1464 Eigen::RowVector3d e1 = V.row(F(i, 2)) - V.row(F(i, 1));
1465 Eigen::RowVector3d e2 = V.row(F(i, 0)) - V.row(F(i, 2));
1466 Eigen::RowVector3d e3 = V.row(F(i, 3)) - V.row(F(i, 0));
1467 Eigen::RowVector3d e4 = V.row(F(i, 3)) - V.row(F(i, 1));
1468 Eigen::RowVector3d e5 = V.row(F(i, 3)) - V.row(F(i, 2));
1469
1470 double l0 = e0.norm();
1471 double l1 = e1.norm();
1472 double l2 = e2.norm();
1473 double l3 = e3.norm();
1474 double l4 = e4.norm();
1475 double l5 = e5.norm();
1476
1477 double J = std::abs((e0.cross(e3)).dot(e2));
1478
1479 double a1 = l0 * l2 * l3;
1480 double a2 = l0 * l1 * l4;
1481 double a3 = l1 * l2 * l5;
1482 double a4 = l3 * l4 * l5;
1483
1484 double a = std::max({a1, a2, a3, a4, J});
1485 quality(i) = J * sqrt(2) / a;
1486 }
1487 }
1488 }
1489
1490} // 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:660
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:440
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