PolyFEM
Loading...
Searching...
No Matches
AdjointTools.cpp
Go to the documentation of this file.
1#include "AdjointTools.hpp"
2
3#include <polyfem/varforms/diff/DifferentiableVarForm.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::varform::DifferentiableVarForm &varform)
71 {
72 if (varform.get_args()["time"]["integrator"].is_string())
73 return 1;
74 if (varform.get_args()["time"]["integrator"]["type"] == "ImplicitEuler")
75 return 1;
76 if (varform.get_args()["time"]["integrator"]["type"] == "BDF")
77 return varform.get_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 varform::DifferentiableVarForm &varform,
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 = varform.primary_space().basis_list();
216 const auto &gbases = varform.primary_space().geometry_basis_list();
217
218 const int dim = varform.get_mesh().dimension();
219 const int actual_dim = varform.get_problem().is_scalar() ? 1 : dim;
220 const int n_elements = int(bases.size());
221 const double t0 = varform.get_problem().is_time_dependent() ? varform.get_args()["time"]["t0"].get<double>() : 0.0;
222 const double dt = varform.get_problem().is_time_dependent() ? varform.get_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(varform.get_mesh().get_body_id(e)) == interested_ids.end())
241 continue;
242
243 assembler::ElementAssemblyValues &vals = local_storage.vals;
244 varform.assembly_cache().compute(e, varform.get_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(varform.primary_assembler().parameters(), e, params.t, quadrature.points, vals.val);
251
252 params.elem = e;
253 params.body_id = varform.get_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(varform.boundary_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 = varform.boundary_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(varform.get_mesh().get_boundary_id(global_primitive_id)) == interested_ids.end())
287 continue;
288
289 utils::BoundarySampler::boundary_quadrature(lb, varform.n_boundary_samples(), varform.get_mesh(), i, false, uv, points, normal, weights);
290
291 assembler::ElementAssemblyValues &vals = local_storage.vals;
292 vals.compute(e, varform.get_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(varform.primary_assembler().parameters(), e, params.t, points, vals.val);
296
297 params.elem = e;
298 params.body_id = varform.get_mesh().get_body_id(e);
299 params.boundary_id = varform.get_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(varform.primary_space().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(varform.primary_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 = varform.get_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 varform::DifferentiableVarForm &varform,
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 = varform.primary_space().geometry_basis_list();
352 const auto &bases = varform.primary_space().basis_list();
353 const int dim = varform.get_mesh().dimension();
354 const int actual_dim = varform.get_problem().is_scalar() ? 1 : dim;
355 const double t0 = varform.get_problem().is_time_dependent() ? varform.get_args()["time"]["t0"].get<double>() : 0.0;
356 const double dt = varform.get_problem().is_time_dependent() ? varform.get_args()["time"]["dt"].get<double>() : 0.0;
357
358 const int n_elements = int(bases.size());
359 term.setZero(varform.primary_space().geometry->n_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(varform.get_mesh().get_body_id(e)) == interested_ids.end())
377 continue;
378
379 assembler::ElementAssemblyValues &vals = local_storage.vals;
380 varform.assembly_cache().compute(e, varform.get_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, varform.get_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(varform.primary_assembler().parameters(), e, params.t, quadrature.points, vals.val);
390
391 params.elem = e;
392 params.body_id = varform.get_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(varform.boundary_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 = varform.boundary_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(varform.get_mesh().get_boundary_id(global_primitive_id)) == interested_ids.end())
455 continue;
456
457 utils::BoundarySampler::boundary_quadrature(lb, varform.n_boundary_samples(), varform.get_mesh(), i, false, uv, points, normal, weights);
458
459 assembler::ElementAssemblyValues &vals = local_storage.vals;
460 io::Evaluator::interpolate_at_local_vals(varform.get_mesh(), varform.get_problem().is_scalar(), bases, gbases, e, points, solution, u, grad_u);
461 // io::Evaluator::interpolate_at_local_vals(varform.get_mesh(), varform.get_problem().is_scalar(), gbases, gbases, e, points, global_positions, x, grad_x);
462
463 vals.compute(e, varform.get_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(varform.primary_assembler().parameters(), e, params.t, points, vals.val);
470
471 params.elem = e;
472 params.body_id = varform.get_mesh().get_body_id(e);
473 params.boundary_id = varform.get_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), varform.get_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 (varform.get_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)(varform.primitive_to_node(), Eigen::all));
574 }
575
576 void AdjointTools::dJ_shape_static_adjoint_term(
577 const varform::DifferentiableVarForm &varform,
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(varform.primary_space().geometry->n_bases * varform.get_mesh().dimension());
586 Eigen::MatrixXd adjoint_zeroed = adjoint;
587 adjoint_zeroed(varform.boundary_state().boundary_nodes, Eigen::all).setZero();
588
589 // if (j.depend_on_u() || j.depend_on_gradu())
590 {
591 ElasticForceDerivative::force_shape_derivative(*varform.solve_data()->elastic_form, 0, varform.primary_space().geometry->n_bases, sol, sol, adjoint_zeroed, elasticity_term);
592 if (varform.solve_data()->body_form)
593 BodyForceDerivative::force_shape_derivative(*varform.solve_data()->body_form, varform.primary_space().geometry->n_bases, 0, sol, adjoint_zeroed, rhs_term);
594 else
595 rhs_term.setZero(one_form.size());
596
597 if (varform.solve_data()->pressure_form)
598 {
599 PressureForceDerivative::force_shape_derivative(*varform.solve_data()->pressure_form, varform.primary_space().geometry->n_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 (varform.is_contact_enabled())
606 {
607 if (const auto barrier_contact = dynamic_cast<const BarrierContactForm *>(varform.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 *>(varform.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 (varform.is_adhesion_enabled())
622 {
623 NormalAdhesionForceDerivative::force_shape_derivative(*varform.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, varform.get_mesh().dimension())(varform.primitive_to_node(), Eigen::all));
634 }
635
636 void AdjointTools::dJ_shape_homogenization_adjoint_term(
637 const varform::DifferentiableVarForm &varform,
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>(varform.solve_data()->nl_problem);
646 assert(homo_problem);
647
648 const int dim = varform.get_mesh().dimension();
649 one_form.setZero(varform.primary_space().geometry->n_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(varform.primary_space().n_bases, varform.primary_space().mesh_nodes, affine_adjoint);
653
654 ElasticForceDerivative::force_shape_derivative(*varform.solve_data()->elastic_form, 0, varform.primary_space().geometry->n_bases, sol, sol, full_adjoint, elasticity_term);
655
656 if (varform.solve_data()->contact_form)
657 {
658 if (const auto barrier_contact = dynamic_cast<const BarrierContactForm *>(varform.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 *>(varform.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 (varform.is_adhesion_enabled())
673 {
674 NormalAdhesionForceDerivative::force_shape_derivative(*varform.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)(varform.primitive_to_node(), Eigen::all));
689 }
690
691 void AdjointTools::dJ_periodic_shape_adjoint_term(
692 const varform::DifferentiableVarForm &varform,
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>(varform.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(varform.primary_space().n_bases, varform.primary_space().mesh_nodes, affine_adjoint);
709
710 const int dim = varform.get_mesh().dimension();
711
712 dJ_shape_homogenization_adjoint_term(varform, 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)(varform.primitive_to_node(), Eigen::all));
720
721 one_form = periodic_mesh_map.apply_jacobian(one_form, periodic_mesh_representation);
722
723 if (varform.solve_data()->periodic_contact_form)
724 {
725 Eigen::VectorXd contact_term;
726 PeriodicContactForceDerivative::force_shape_derivative(*varform.solve_data()->periodic_contact_form, varform, diff_cache, periodic_mesh_map, periodic_mesh_representation, varform.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 varform::DifferentiableVarForm &varform,
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 = varform.get_args()["time"]["t0"];
740 const double dt = varform.get_args()["time"]["dt"];
741 const int time_steps = varform.get_args()["time"]["time_steps"];
742 const int bdf_order = get_bdf_order(varform);
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(varform.primary_space().geometry->n_bases * varform.get_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(varform.boundary_state().boundary_nodes).setZero();
760 cur_nu(varform.boundary_state().boundary_nodes).setZero();
761
762 {
763 InertiaForceDerivative::force_shape_derivative(*varform.solve_data()->inertia_form, varform.get_mesh().is_volume(), varform.primary_space().geometry->n_bases, t, varform.primary_space().basis_list(), varform.primary_space().geometry_basis_list(), varform.mass_assembler(), varform.mass_assembly_cache(), velocity, cur_nu, mass_term);
764 ElasticForceDerivative::force_shape_derivative(*varform.solve_data()->elastic_form, t, varform.primary_space().geometry->n_bases, diff_cache.u(i), diff_cache.u(i), cur_p, elasticity_term);
765 BodyForceDerivative::force_shape_derivative(*varform.solve_data()->body_form, varform.primary_space().geometry->n_bases, t, diff_cache.u(i - 1), cur_p, rhs_term);
766 PressureForceDerivative::force_shape_derivative(*varform.solve_data()->pressure_form, varform.primary_space().geometry->n_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 (varform.solve_data()->damping_form)
770 ElasticForceDerivative::force_shape_derivative(*varform.solve_data()->damping_form, t, varform.primary_space().geometry->n_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 (varform.is_contact_enabled())
775 {
776 if (const auto barrier_contact = dynamic_cast<const BarrierContactForm *>(varform.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 *>(varform.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 (varform.solve_data()->friction_form)
791 {
792 FrictionForceDerivative::force_shape_derivative(*varform.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 (varform.is_adhesion_enabled())
800 {
801 NormalAdhesionForceDerivative::force_shape_derivative(*varform.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 (varform.solve_data()->tangential_adhesion_form)
810 {
811 TangentialAdhesionForceDerivative::force_shape_derivative(*varform.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(varform.boundary_state().boundary_nodes).setZero();
834 InertiaForceDerivative::force_shape_derivative(*varform.solve_data()->inertia_form, varform.get_mesh().is_volume(), varform.primary_space().geometry->n_bases, t0, varform.primary_space().basis_list(), varform.primary_space().geometry_basis_list(), varform.mass_assembler(), varform.mass_assembly_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, varform.get_mesh().dimension())(varform.primitive_to_node(), Eigen::all));
839 }
840
841 void AdjointTools::dJ_material_static_adjoint_term(
842 const varform::DifferentiableVarForm &varform,
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(varform.boundary_state().boundary_nodes, Eigen::all).setZero();
849 ElasticForceDerivative::force_material_derivative(*varform.solve_data()->elastic_form, 0, sol, sol, adjoint_zeroed, one_form);
850 }
851
852 void AdjointTools::dJ_material_transient_adjoint_term(
853 const varform::DifferentiableVarForm &varform,
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 = varform.get_args()["time"]["t0"];
860 const double dt = varform.get_args()["time"]["dt"];
861 const int time_steps = varform.get_args()["time"]["time_steps"];
862 const int bdf_order = get_bdf_order(varform);
863
864 one_form.setZero(varform.primary_space().basis_list().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(varform.boundary_state().boundary_nodes).setZero();
879
880 ElasticForceDerivative::force_material_derivative(*varform.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 varform::DifferentiableVarForm &varform,
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 = varform.get_args()["time"]["dt"];
897 const double mu = varform.solve_data()->friction_form->mu();
898 const int time_steps = varform.get_args()["time"]["time_steps"];
899 const int dim = varform.get_mesh().dimension();
900 const int bdf_order = get_bdf_order(varform);
901
902 one_form.setZero(1);
903
904 std::shared_ptr<time_integrator::ImplicitTimeIntegrator> time_integrator =
905 time_integrator::ImplicitTimeIntegrator::construct_time_integrator(varform.get_args()["time"]["integrator"]);
906 {
907 Eigen::MatrixXd solution, velocity, acceleration;
908
909 const varform::InitialConditionOverride *ic_override =
910 diff_cache.initial_condition_override ? &*diff_cache.initial_condition_override : nullptr;
911
912 solution = diff_cache.u(0);
913 varform.initial_velocity(velocity, ic_override);
914 varform.initial_acceleration(acceleration, ic_override);
915 const double dt = varform.get_args()["time"]["dt"];
916 time_integrator->init(solution, velocity, acceleration, dt);
917 }
918
919 for (int t = 1; t <= time_steps; ++t)
920 {
921 const int real_order = std::min(bdf_order, t);
922 double beta = time_integrator::BDF::betas(real_order - 1);
923
924 const Eigen::MatrixXd surface_solution_prev = varform.collision_mesh().vertices(utils::unflatten(diff_cache.u(t - 1), dim));
925 // const Eigen::MatrixXd surface_solution = varform.collision_mesh().vertices(utils::unflatten(diff_cache.u(t), dim));
926
927 const Eigen::MatrixXd surface_velocities = varform.collision_mesh().map_displacements(utils::unflatten(time_integrator->compute_velocity(diff_cache.u(t)), varform.collision_mesh().dim()));
928 time_integrator->update_quantities(diff_cache.u(t));
929
930 if (const auto barrier_contact = dynamic_cast<const BarrierContactForm *>(varform.solve_data()->contact_form.get()))
931 {
932 ipc::BarrierPotential bp = barrier_contact->barrier_potential();
933 bp.set_stiffness(barrier_contact->barrier_stiffness());
934 Eigen::MatrixXd force = varform.collision_mesh().to_full_dof(
935 -varform.solve_data()->friction_form->friction_potential().force(
936 diff_cache.friction_collision_set(t),
937 varform.collision_mesh(),
938 varform.collision_mesh().rest_positions(),
939 /*lagged_displacements=*/surface_solution_prev,
940 surface_velocities,
941 bp,
942 0., true));
943
944 Eigen::VectorXd cur_p = adjoint_p.col(t);
945 cur_p(varform.boundary_state().boundary_nodes).setZero();
946
947 one_form(0) += dot(cur_p, force) * beta * dt;
948 }
949 }
950 }
951
952 void AdjointTools::dJ_damping_transient_adjoint_term(
953 const varform::DifferentiableVarForm &varform,
954 const DiffCache &diff_cache,
955 const Eigen::MatrixXd &adjoint_nu,
956 const Eigen::MatrixXd &adjoint_p,
957 Eigen::VectorXd &one_form)
958 {
959 const double t0 = varform.get_args()["time"]["t0"];
960 const double dt = varform.get_args()["time"]["dt"];
961 const int time_steps = varform.get_args()["time"]["time_steps"];
962 const int bdf_order = get_bdf_order(varform);
963
964 one_form.setZero(2);
965
966 auto storage = utils::create_thread_storage(LocalThreadVecStorage(one_form.size()));
967
968 utils::maybe_parallel_for(time_steps, [&](int start, int end, int thread_id) {
969 LocalThreadVecStorage &local_storage = utils::get_local_thread_storage(storage, thread_id);
970 Eigen::VectorXd damping_term;
971 for (int t_aux = start; t_aux < end; ++t_aux)
972 {
973 const int t = time_steps - t_aux;
974 const int real_order = std::min(bdf_order, t);
975 const double beta = time_integrator::BDF::betas(real_order - 1);
976
977 Eigen::VectorXd cur_p = adjoint_p.col(t);
978 cur_p(varform.boundary_state().boundary_nodes).setZero();
979
980 ElasticForceDerivative::force_material_derivative(*varform.solve_data()->damping_form, t * dt + t0, diff_cache.u(t), diff_cache.u(t - 1), -cur_p, damping_term);
981 local_storage.vec += (beta * dt) * damping_term;
982 }
983 });
984
985 for (const LocalThreadVecStorage &local_storage : storage)
986 one_form += local_storage.vec;
987 }
988
989 void AdjointTools::dJ_initial_condition_adjoint_term(
990 const varform::DifferentiableVarForm &varform,
991 const Eigen::MatrixXd &adjoint_nu,
992 const Eigen::MatrixXd &adjoint_p,
993 Eigen::VectorXd &one_form)
994 {
995 const int ndof = varform.primary_space().ndof();
996 one_form.setZero(ndof * 2); // half for initial solution, half for initial velocity
997
998 // \partial_q \hat{J}^0 - p_0^T \partial_q g^v - \mu_0^T \partial_q g^u
999 one_form.segment(0, ndof) = -adjoint_nu.col(0); // adjoint_nu[0] actually stores adjoint_mu[0]
1000 one_form.segment(ndof, ndof) = -adjoint_p.col(0);
1001
1002 for (int b : varform.boundary_state().boundary_nodes)
1003 {
1004 one_form(b) = 0;
1005 one_form(ndof + b) = 0;
1006 }
1007 }
1008
1009 void AdjointTools::dJ_dirichlet_static_adjoint_term(
1010 const varform::DifferentiableVarForm &varform,
1011 const DiffCache &diff_cache,
1012 const Eigen::MatrixXd &adjoint,
1013 Eigen::VectorXd &one_form)
1014 {
1015 const int n_dirichlet_dof = varform.boundary_state().boundary_nodes.size();
1016 StiffnessMatrix gradd_h = diff_cache.gradu_h(0);
1017 std::set<int> boundary_nodes_set(varform.boundary_state().boundary_nodes.begin(), varform.boundary_state().boundary_nodes.end());
1018 gradd_h.prune([&boundary_nodes_set](const Eigen::Index &row, const Eigen::Index &col, const FullNLProblem::Scalar &value) {
1019 if (row != col)
1020 return value;
1021 if (boundary_nodes_set.find(row) == boundary_nodes_set.end())
1022 return value;
1023 return 0.0;
1024 });
1025 one_form.setZero(varform.primary_space().ndof());
1026 one_form(varform.boundary_state().boundary_nodes) -= adjoint(varform.boundary_state().boundary_nodes, 0);
1027 one_form(varform.boundary_state().boundary_nodes) -= (gradd_h.transpose() * adjoint.col(0))(varform.boundary_state().boundary_nodes);
1028 one_form = utils::flatten(utils::unflatten(one_form, varform.get_mesh().dimension())(varform.primitive_to_node(), Eigen::all));
1029 }
1030
1031 void AdjointTools::dJ_dirichlet_transient_adjoint_term(
1032 const varform::DifferentiableVarForm &varform,
1033 const Eigen::MatrixXd &adjoint_nu,
1034 const Eigen::MatrixXd &adjoint_p,
1035 Eigen::VectorXd &one_form)
1036 {
1037 const double dt = varform.get_args()["time"]["dt"];
1038 const int time_steps = varform.get_args()["time"]["time_steps"];
1039 const int bdf_order = get_bdf_order(varform);
1040 const int n_dirichlet_dof = varform.boundary_state().boundary_nodes.size();
1041
1042 // Map dirichlet gradient on each node to dirichlet gradient on boundary ids
1043
1044 one_form.setZero(time_steps * n_dirichlet_dof);
1045 for (int i = 1; i <= time_steps; ++i)
1046 {
1047 const int real_order = std::min(bdf_order, i);
1048 const double beta_dt = time_integrator::BDF::betas(real_order - 1) * dt;
1049
1050 one_form.segment((i - 1) * n_dirichlet_dof, n_dirichlet_dof) = -(1. / beta_dt) * adjoint_p(varform.boundary_state().boundary_nodes, i);
1051 }
1052 }
1053
1054 void AdjointTools::dJ_pressure_static_adjoint_term(
1055 const varform::DifferentiableVarForm &varform,
1056 const std::vector<int> &boundary_ids,
1057 const Eigen::MatrixXd &sol,
1058 const Eigen::MatrixXd &adjoint,
1059 Eigen::VectorXd &one_form)
1060 {
1061 const int n_pressure_dof = boundary_ids.size();
1062
1063 one_form.setZero(n_pressure_dof);
1064
1065 for (int i = 0; i < boundary_ids.size(); ++i)
1066 {
1067 double pressure_term = PressureForceDerivative::force_pressure_derivative(
1068 *varform.solve_data()->pressure_form,
1069 varform.primary_space().geometry->n_bases,
1070 0,
1071 boundary_ids[i],
1072 sol,
1073 adjoint);
1074 one_form(i) = pressure_term;
1075 }
1076 }
1077
1078 void AdjointTools::dJ_pressure_transient_adjoint_term(
1079 const varform::DifferentiableVarForm &varform,
1080 const DiffCache &diff_cache,
1081 const std::vector<int> &boundary_ids,
1082 const Eigen::MatrixXd &adjoint_nu,
1083 const Eigen::MatrixXd &adjoint_p,
1084 Eigen::VectorXd &one_form)
1085 {
1086 const double t0 = varform.get_args()["time"]["t0"];
1087 const double dt = varform.get_args()["time"]["dt"];
1088 const int time_steps = varform.get_args()["time"]["time_steps"];
1089 const int bdf_order = get_bdf_order(varform);
1090
1091 const int n_pressure_dof = boundary_ids.size();
1092
1093 one_form.setZero(time_steps * n_pressure_dof);
1094 Eigen::VectorXd cur_p, cur_nu;
1095 for (int i = time_steps; i > 0; --i)
1096 {
1097 const int real_order = std::min(bdf_order, i);
1098 double beta = time_integrator::BDF::betas(real_order - 1);
1099 double beta_dt = beta * dt;
1100 const double t = i * dt + t0;
1101
1102 cur_p = adjoint_p.col(i);
1103 cur_nu = adjoint_nu.col(i);
1104 cur_p(varform.boundary_state().boundary_nodes).setZero();
1105 cur_nu(varform.boundary_state().boundary_nodes).setZero();
1106
1107 for (int b = 0; b < boundary_ids.size(); ++b)
1108 {
1109 double pressure_term = PressureForceDerivative::force_pressure_derivative(
1110 *varform.solve_data()->pressure_form,
1111 varform.primary_space().geometry->n_bases,
1112 t,
1113 boundary_ids[b],
1114 diff_cache.u(i),
1115 cur_p);
1116 one_form((i - 1) * n_pressure_dof + b) = -beta_dt * pressure_term;
1117 }
1118 }
1119 }
1120
1121 void AdjointTools::dJ_du_step(
1122 const varform::DifferentiableVarForm &varform,
1123 const IntegrableFunctional &j,
1124 const Eigen::MatrixXd &solution,
1125 const std::set<int> &interested_ids,
1126 const SpatialIntegralType spatial_integral_type,
1127 const int cur_step,
1128 Eigen::VectorXd &term)
1129 {
1130 const auto &bases = varform.primary_space().basis_list();
1131 const auto &gbases = varform.primary_space().geometry_basis_list();
1132
1133 const int dim = varform.get_mesh().dimension();
1134 const int actual_dim = varform.get_problem().is_scalar() ? 1 : dim;
1135 const int n_elements = int(bases.size());
1136 const double t0 = varform.get_problem().is_time_dependent() ? varform.get_args()["time"]["t0"].get<double>() : 0.0;
1137 const double dt = varform.get_problem().is_time_dependent() ? varform.get_args()["time"]["dt"].get<double>() : 0.0;
1138
1139 term = Eigen::MatrixXd::Zero(varform.primary_space().n_bases * actual_dim, 1);
1140
1141 if (!j.depend_on_u() && !j.depend_on_gradu() && !j.depend_on_gradu_local())
1142 return;
1143
1144 if (spatial_integral_type == SpatialIntegralType::Volume)
1145 {
1146 auto storage = utils::create_thread_storage(LocalThreadVecStorage(term.size()));
1147 utils::maybe_parallel_for(n_elements, [&](int start, int end, int thread_id) {
1148 LocalThreadVecStorage &local_storage = utils::get_local_thread_storage(storage, thread_id);
1149
1150 Eigen::MatrixXd u, grad_u;
1151 Eigen::MatrixXd lambda, mu;
1152 Eigen::MatrixXd dj_du, dj_dgradu, dj_dgradx;
1153
1154 IntegrableFunctional::ParameterType params;
1155 params.t = dt * cur_step + t0;
1156 params.step = cur_step;
1157
1158 for (int e = start; e < end; ++e)
1159 {
1160 if (interested_ids.size() != 0 && interested_ids.find(varform.get_mesh().get_body_id(e)) == interested_ids.end())
1161 continue;
1162
1163 assembler::ElementAssemblyValues &vals = local_storage.vals;
1164 varform.assembly_cache().compute(e, varform.get_mesh().is_volume(), bases[e], gbases[e], vals);
1165
1166 const quadrature::Quadrature &quadrature = vals.quadrature;
1167 local_storage.da = vals.det.array() * quadrature.weights.array();
1168
1169 const Eigen::MatrixXd lame_params = extract_lame_params(varform.primary_assembler().parameters(), e, params.t, quadrature.points, vals.val);
1170
1171 const int n_loc_bases_ = int(vals.basis_values.size());
1172
1173 io::Evaluator::interpolate_at_local_vals(e, dim, actual_dim, vals, solution, u, grad_u);
1174
1175 params.elem = e;
1176 params.body_id = varform.get_mesh().get_body_id(e);
1177
1178 dj_dgradu.resize(0, 0);
1179 if (j.depend_on_gradu())
1180 {
1181 j.dj_dgradu(lame_params, quadrature.points, vals.val, u, grad_u, Eigen::MatrixXd::Zero(0, 0) /*Not used*/, vals, params, dj_dgradu);
1182 for (int q = 0; q < dj_dgradu.rows(); q++)
1183 dj_dgradu.row(q) *= local_storage.da(q);
1184 }
1185
1186 dj_du.resize(0, 0);
1187 if (j.depend_on_u())
1188 {
1189 j.dj_du(lame_params, quadrature.points, vals.val, u, grad_u, Eigen::MatrixXd::Zero(0, 0) /*Not used*/, vals, params, dj_du);
1190 for (int q = 0; q < dj_du.rows(); q++)
1191 dj_du.row(q) *= local_storage.da(q);
1192 }
1193
1194 for (int i = 0; i < n_loc_bases_; ++i)
1195 {
1196 const assembler::AssemblyValues &v = vals.basis_values[i];
1197 assert(v.global.size() == 1);
1198 for (int d = 0; d < actual_dim; d++)
1199 {
1200 double val = 0;
1201
1202 // j = j(..., grad u)
1203 if (j.depend_on_gradu())
1204 {
1205 for (int q = 0; q < local_storage.da.size(); ++q)
1206 val += dot(dj_dgradu.block(q, d * dim, 1, dim), v.grad_t_m.row(q));
1207 }
1208
1209 // j = j(..., u)
1210 if (j.depend_on_u())
1211 {
1212 for (int q = 0; q < local_storage.da.size(); ++q)
1213 val += dj_du(q, d) * v.val(q);
1214 }
1215 local_storage.vec(v.global[0].index * actual_dim + d) += val;
1216 }
1217 }
1218 }
1219 });
1220 for (const LocalThreadVecStorage &local_storage : storage)
1221 term += local_storage.vec;
1222 }
1223 else if (spatial_integral_type == SpatialIntegralType::Surface)
1224 {
1225 auto storage = utils::create_thread_storage(LocalThreadVecStorage(term.size()));
1226 utils::maybe_parallel_for(varform.boundary_state().total_local_boundary.size(), [&](int start, int end, int thread_id) {
1227 LocalThreadVecStorage &local_storage = utils::get_local_thread_storage(storage, thread_id);
1228
1229 Eigen::MatrixXd uv, samples, gtmp;
1230 Eigen::MatrixXd points, normal;
1231 Eigen::VectorXd weights;
1232
1233 Eigen::MatrixXd u, grad_u;
1234 Eigen::MatrixXd lambda, mu;
1235 Eigen::MatrixXd dj_du, dj_dgradu, dj_dgradu_local;
1236
1237 IntegrableFunctional::ParameterType params;
1238 params.t = dt * cur_step + t0;
1239 params.step = cur_step;
1240
1241 for (int lb_id = start; lb_id < end; ++lb_id)
1242 {
1243 const auto &lb = varform.boundary_state().total_local_boundary[lb_id];
1244 const int e = lb.element_id();
1245
1246 for (int i = 0; i < lb.size(); i++)
1247 {
1248 const int global_primitive_id = lb.global_primitive_id(i);
1249 if (interested_ids.size() != 0 && interested_ids.find(varform.get_mesh().get_boundary_id(global_primitive_id)) == interested_ids.end())
1250 continue;
1251
1252 utils::BoundarySampler::boundary_quadrature(lb, varform.n_boundary_samples(), varform.get_mesh(), i, false, uv, points, normal, weights);
1253
1254 assembler::ElementAssemblyValues &vals = local_storage.vals;
1255 vals.compute(e, varform.get_mesh().is_volume(), points, bases[e], gbases[e]);
1256 io::Evaluator::interpolate_at_local_vals(e, dim, actual_dim, vals, solution, u, grad_u);
1257
1258 const Eigen::MatrixXd lame_params = extract_lame_params(varform.primary_assembler().parameters(), e, params.t, points, vals.val);
1259
1260 // normal = normal * vals.jac_it[0]; // assuming linear geometry
1261
1262 const int n_loc_bases_ = int(vals.basis_values.size());
1263
1264 params.elem = e;
1265 params.body_id = varform.get_mesh().get_body_id(e);
1266 params.boundary_id = varform.get_mesh().get_boundary_id(global_primitive_id);
1267
1268 dj_dgradu.resize(0, 0);
1269 if (j.depend_on_gradu())
1270 {
1271 j.dj_dgradu(lame_params, points, vals.val, u, grad_u, normal, vals, params, dj_dgradu);
1272 for (int q = 0; q < dj_dgradu.rows(); q++)
1273 dj_dgradu.row(q) *= weights(q);
1274 }
1275
1276 dj_dgradu_local.resize(0, 0);
1277 if (j.depend_on_gradu_local())
1278 {
1279 j.dj_dgradu_local(lame_params, points, vals.val, u, grad_u, normal, vals, params, dj_dgradu_local);
1280 for (int q = 0; q < dj_dgradu_local.rows(); q++)
1281 dj_dgradu_local.row(q) *= weights(q);
1282 }
1283
1284 dj_du.resize(0, 0);
1285 if (j.depend_on_u())
1286 {
1287 j.dj_du(lame_params, points, vals.val, u, grad_u, normal, vals, params, dj_du);
1288 for (int q = 0; q < dj_du.rows(); q++)
1289 dj_du.row(q) *= weights(q);
1290 }
1291
1292 for (int l = 0; l < lb.size(); ++l)
1293 {
1294 const auto nodes = bases[e].local_nodes_for_primitive(lb.global_primitive_id(l), varform.get_mesh());
1295
1296 for (long n = 0; n < nodes.size(); ++n)
1297 {
1298 const assembler::AssemblyValues &v = vals.basis_values[nodes(n)];
1299 assert(v.global.size() == 1);
1300 for (int d = 0; d < actual_dim; d++)
1301 {
1302 double val = 0;
1303
1304 // j = j(x, grad u)
1305 if (j.depend_on_gradu())
1306 {
1307 for (int q = 0; q < weights.size(); ++q)
1308 val += dot(dj_dgradu.block(q, d * dim, 1, dim), v.grad_t_m.row(q));
1309 }
1310 // j = j(x, grad u)
1311 if (j.depend_on_gradu_local())
1312 {
1313 for (int q = 0; q < weights.size(); ++q)
1314 val += dot(dj_dgradu_local.block(q, d * dim, 1, dim), v.grad.row(q));
1315 }
1316 // j = j(x, u)
1317 if (j.depend_on_u())
1318 {
1319 for (int q = 0; q < weights.size(); ++q)
1320 val += dj_du(q, d) * v.val(q);
1321 }
1322 local_storage.vec(v.global[0].index * actual_dim + d) += val;
1323 }
1324 }
1325 }
1326 }
1327 }
1328 });
1329 for (const LocalThreadVecStorage &local_storage : storage)
1330 term += local_storage.vec;
1331 }
1332 else if (spatial_integral_type == SpatialIntegralType::VertexSum)
1333 {
1334 std::vector<bool> traversed(varform.primary_space().n_bases, false);
1335 IntegrableFunctional::ParameterType params;
1336 params.t = dt * cur_step + t0;
1337 params.step = cur_step;
1338 for (int e = 0; e < bases.size(); e++)
1339 {
1340 const auto &bs = bases[e];
1341 for (int i = 0; i < bs.bases.size(); i++)
1342 {
1343 const auto &b = bs.bases[i];
1344 assert(b.global().size() == 1);
1345 const auto &g = b.global()[0];
1346 if (traversed[g.index])
1347 continue;
1348
1349 const Eigen::MatrixXd lame_params = extract_lame_params(varform.primary_assembler().parameters(), e, params.t, Eigen::MatrixXd::Zero(1, dim) /*Not used*/, g.node);
1350
1351 params.node = g.index;
1352 params.elem = e;
1353 params.body_id = varform.get_mesh().get_body_id(e);
1354 Eigen::MatrixXd val;
1355 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);
1356 term.block(g.index * actual_dim, 0, actual_dim, 1) += val.transpose();
1357 traversed[g.index] = true;
1358 }
1359 }
1360 }
1361 }
1362
1363 Eigen::VectorXd AdjointTools::map_primitive_to_node_order(const varform::DifferentiableVarForm &varform, const Eigen::VectorXd &primitives)
1364 {
1365 int dim = varform.get_mesh().dimension();
1366 assert(primitives.size() == (varform.primary_space().geometry->n_bases * dim));
1367 Eigen::VectorXd nodes(primitives.size());
1368 auto map = varform.primitive_to_node();
1369 for (int v = 0; v < varform.primary_space().geometry->n_bases; ++v)
1370 nodes.segment(map[v] * dim, dim) = primitives.segment(v * dim, dim);
1371 return nodes;
1372 }
1373
1374 Eigen::VectorXd AdjointTools::map_node_to_primitive_order(const varform::DifferentiableVarForm &varform, const Eigen::VectorXd &nodes)
1375 {
1376 int dim = varform.get_mesh().dimension();
1377 assert(nodes.size() == (varform.primary_space().geometry->n_bases * dim));
1378 Eigen::VectorXd primitives(nodes.size());
1379 auto map = varform.node_to_primitive();
1380 for (int v = 0; v < varform.primary_space().geometry->n_bases; ++v)
1381 primitives.segment(map[v] * dim, dim) = nodes.segment(v * dim, dim);
1382 return primitives;
1383 }
1384
1385 Eigen::MatrixXd AdjointTools::edge_normal_gradient(const Eigen::MatrixXd &V)
1386 {
1387 DiffScalarBase::setVariableCount(4);
1388 Eigen::Matrix<Diff, 4, 1> full_diff(4, 1);
1389 for (int i = 0; i < 2; i++)
1390 for (int j = 0; j < 2; j++)
1391 full_diff(i * 2 + j) = Diff(i * 2 + j, V(i, j));
1392 auto reduced_diff = edge_normal(full_diff);
1393
1394 Eigen::MatrixXd grad(2, 4);
1395 for (int i = 0; i < 2; ++i)
1396 grad.row(i) = reduced_diff[i].getGradient();
1397
1398 return grad;
1399 }
1400
1401 Eigen::MatrixXd AdjointTools::face_normal_gradient(const Eigen::MatrixXd &V)
1402 {
1403 DiffScalarBase::setVariableCount(9);
1404 Eigen::Matrix<Diff, 9, 1> full_diff(9, 1);
1405 for (int i = 0; i < 3; i++)
1406 for (int j = 0; j < 3; j++)
1407 full_diff(i * 3 + j) = Diff(i * 3 + j, V(i, j));
1408 auto reduced_diff = face_normal(full_diff);
1409
1410 Eigen::MatrixXd grad(3, 9);
1411 for (int i = 0; i < 3; ++i)
1412 grad.row(i) = reduced_diff[i].getGradient();
1413
1414 return grad;
1415 }
1416
1417 Eigen::MatrixXd AdjointTools::edge_velocity_divergence(const Eigen::MatrixXd &V)
1418 {
1419 return line_length_grad(V) / line_length<double>(V);
1420 }
1421
1422 Eigen::MatrixXd AdjointTools::face_velocity_divergence(const Eigen::MatrixXd &V)
1423 {
1424 return triangle_area_grad(V) / triangle_area<double>(V);
1425 }
1426
1427 void AdjointTools::scaled_jacobian(const Eigen::MatrixXd &V, const Eigen::MatrixXi &F, Eigen::VectorXd &quality)
1428 {
1429 const int dim = F.cols() - 1;
1430
1431 quality.setZero(F.rows());
1432 if (dim == 2)
1433 {
1434 for (int i = 0; i < F.rows(); i++)
1435 {
1436 Eigen::RowVector3d e0;
1437 e0(2) = 0;
1438 e0.head(2) = V.row(F(i, 2)) - V.row(F(i, 1));
1439 Eigen::RowVector3d e1;
1440 e1(2) = 0;
1441 e1.head(2) = V.row(F(i, 0)) - V.row(F(i, 2));
1442 Eigen::RowVector3d e2;
1443 e2(2) = 0;
1444 e2.head(2) = V.row(F(i, 1)) - V.row(F(i, 0));
1445
1446 double l0 = e0.norm();
1447 double l1 = e1.norm();
1448 double l2 = e2.norm();
1449
1450 double A = 0.5 * (e0.cross(e1)).norm();
1451 double Lmax = std::max(l0 * l1, std::max(l1 * l2, l0 * l2));
1452
1453 quality(i) = 2 * A * (2 / sqrt(3)) / Lmax;
1454 }
1455 }
1456 else
1457 {
1458 for (int i = 0; i < F.rows(); i++)
1459 {
1460 Eigen::RowVector3d e0 = V.row(F(i, 1)) - V.row(F(i, 0));
1461 Eigen::RowVector3d e1 = V.row(F(i, 2)) - V.row(F(i, 1));
1462 Eigen::RowVector3d e2 = V.row(F(i, 0)) - V.row(F(i, 2));
1463 Eigen::RowVector3d e3 = V.row(F(i, 3)) - V.row(F(i, 0));
1464 Eigen::RowVector3d e4 = V.row(F(i, 3)) - V.row(F(i, 1));
1465 Eigen::RowVector3d e5 = V.row(F(i, 3)) - V.row(F(i, 2));
1466
1467 double l0 = e0.norm();
1468 double l1 = e1.norm();
1469 double l2 = e2.norm();
1470 double l3 = e3.norm();
1471 double l4 = e4.norm();
1472 double l5 = e5.norm();
1473
1474 double J = std::abs((e0.cross(e3)).dot(e2));
1475
1476 double a1 = l0 * l2 * l3;
1477 double a2 = l0 * l1 * l4;
1478 double a3 = l1 * l2 * l5;
1479 double a4 = l3 * l4 * l5;
1480
1481 double a = std::max({a1, a2, a3, a4, J});
1482 quality(i) = J * sqrt(2) / a;
1483 }
1484 }
1485 }
1486
1487} // namespace polyfem::solver
V
int V
Definition AdjointNLProblem.cpp:48
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:76
val
double val
Definition Assembler.cpp:90
da
QuadratureVector da
Definition Assembler.cpp:27
vals
ElementAssemblyValues vals
Definition Assembler.cpp:26
gvals
assembler::ElementAssemblyValues gvals
Definition BodyForceDerivative.cpp:26
BodyForceDerivative.hpp
BodyForm.hpp
quadrature
Quadrature quadrature
Definition MassMatrixAssembler.cpp:30
J
double J
Definition NavierStokesFSI.cpp:60
F
Eigen::MatrixXd F
Definition NavierStokesFSI.cpp:58
polyfem::DiffCache
Storage for additional data required by differntial code.
Definition DiffCache.hpp:22
polyfem::DiffCache::collision_set
const ipc::NormalCollisions & collision_set(int step) const
Definition DiffCache.hpp:94
polyfem::DiffCache::disp_grad
Eigen::MatrixXd disp_grad(int step=0) const
Definition DiffCache.hpp:55
polyfem::DiffCache::initial_condition_override
std::optional< varform::InitialConditionOverride > initial_condition_override
Initial-condition override storage for initial condition optimization.
Definition DiffCache.hpp:25
polyfem::DiffCache::v
Eigen::VectorXd v(int step) const
Definition DiffCache.hpp:70
polyfem::DiffCache::friction_collision_set
const ipc::TangentialCollisions & friction_collision_set(int step) const
Definition DiffCache.hpp:108
polyfem::DiffCache::basis_nodes_to_gbasis_nodes
const StiffnessMatrix & basis_nodes_to_gbasis_nodes() const
Definition DiffCache.cpp:279
polyfem::DiffCache::smooth_collision_set
const ipc::SmoothCollisions & smooth_collision_set(int step) const
Definition DiffCache.hpp:101
polyfem::DiffCache::u
Eigen::VectorXd u(int step) const
Definition DiffCache.hpp:63
polyfem::DiffCache::gradu_h
const StiffnessMatrix & gradu_h(int step) const
Definition DiffCache.hpp:85
polyfem::DiffCache::normal_adhesion_collision_set
const ipc::NormalCollisions & normal_adhesion_collision_set(int step) const
Definition DiffCache.hpp:115
polyfem::DiffCache::tangential_adhesion_collision_set
const ipc::TangentialCollisions & tangential_adhesion_collision_set(int step) const
Definition DiffCache.hpp:122
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::assembler::Assembler::parameters
virtual std::map< std::string, ParamFunc > parameters() const =0
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::assembler::Problem::is_scalar
virtual bool is_scalar() const =0
polyfem::assembler::Problem::is_time_dependent
virtual bool is_time_dependent() const
Definition Problem.hpp:62
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:689
polyfem::mesh::Mesh::get_body_id
virtual int get_body_id(const int primitive) const
Get the volume selection of an element (cell in 3d, face in 2d)
Definition Mesh.hpp:525
polyfem::mesh::Mesh::is_volume
virtual bool is_volume() const =0
checks if mesh is volume
polyfem::mesh::Mesh::dimension
int dimension() const
utily for dimension
Definition Mesh.hpp:164
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:184
polyfem::solver::SolveData::inertia_form
std::shared_ptr< solver::InertiaForm > inertia_form
Definition SolveData.hpp:185
polyfem::solver::SolveData::periodic_contact_form
std::shared_ptr< solver::PeriodicContactForm > periodic_contact_form
Definition SolveData.hpp:190
polyfem::solver::SolveData::pressure_form
std::shared_ptr< solver::PressureForm > pressure_form
Definition SolveData.hpp:186
polyfem::solver::SolveData::body_form
std::shared_ptr< solver::BodyForm > body_form
Definition SolveData.hpp:180
polyfem::solver::SolveData::nl_problem
std::shared_ptr< solver::NLProblem > nl_problem
Definition SolveData.hpp:176
polyfem::solver::SolveData::normal_adhesion_form
std::shared_ptr< solver::NormalAdhesionForm > normal_adhesion_form
Definition SolveData.hpp:187
polyfem::solver::SolveData::contact_form
std::shared_ptr< solver::ContactForm > contact_form
Definition SolveData.hpp:181
polyfem::solver::SolveData::damping_form
std::shared_ptr< solver::ElasticForm > damping_form
Definition SolveData.hpp:182
polyfem::solver::SolveData::elastic_form
std::shared_ptr< solver::ElasticForm > elastic_form
Definition SolveData.hpp:183
polyfem::solver::SolveData::tangential_adhesion_form
std::shared_ptr< solver::TangentialAdhesionForm > tangential_adhesion_form
Definition SolveData.hpp:188
polyfem::varform::DifferentiableVarForm
Optimization-facing interface implemented by differentiated VarForm adapters.
Definition DifferentiableVarForm.hpp:23
polyfem::varform::DifferentiableVarForm::assembly_cache
virtual const assembler::AssemblyValsCache & assembly_cache() const =0
polyfem::varform::DifferentiableVarForm::primary_assembler
virtual const assembler::Assembler & primary_assembler() const =0
polyfem::varform::DifferentiableVarForm::boundary_state
virtual const VarFormBoundaryState & boundary_state() const =0
polyfem::varform::DifferentiableVarForm::get_mesh
virtual const mesh::Mesh & get_mesh() const =0
polyfem::varform::DifferentiableVarForm::solve_data
virtual solver::SolveData * solve_data()=0
polyfem::varform::DifferentiableVarForm::get_problem
virtual assembler::Problem & get_problem()=0
polyfem::varform::DifferentiableVarForm::is_contact_enabled
virtual bool is_contact_enabled() const =0
polyfem::varform::DifferentiableVarForm::primary_space
virtual const FESpace & primary_space() const =0
polyfem::varform::DifferentiableVarForm::mass_assembler
virtual const assembler::Mass & mass_assembler() const =0
polyfem::varform::DifferentiableVarForm::mass_assembly_cache
virtual const assembler::AssemblyValsCache & mass_assembly_cache() const =0
polyfem::varform::DifferentiableVarForm::collision_mesh
virtual const ipc::CollisionMesh & collision_mesh() const
Definition DifferentiableVarForm.cpp:11
polyfem::varform::DifferentiableVarForm::initial_acceleration
virtual void initial_acceleration(Eigen::MatrixXd &acceleration, const InitialConditionOverride *override=nullptr) const
Definition DifferentiableVarForm.cpp:31
polyfem::varform::DifferentiableVarForm::initial_velocity
virtual void initial_velocity(Eigen::MatrixXd &velocity, const InitialConditionOverride *override=nullptr) const
Definition DifferentiableVarForm.cpp:26
polyfem::varform::FESpace::geometry_basis_list
const std::vector< basis::ElementBases > & geometry_basis_list() const
Definition FESpace.hpp:115
polyfem::varform::FESpace::geometry
std::shared_ptr< GeometryMapping > geometry
Geometric mapping used to integrate this FE space.
Definition FESpace.hpp:89
polyfem::varform::FESpace::n_bases
int n_bases
Number of globally indexed scalar basis functions in the space.
Definition FESpace.hpp:65
polyfem::varform::FESpace::basis_list
const std::vector< basis::ElementBases > & basis_list() const
Definition FESpace.hpp:109
polyfem::varform::FESpace::mesh_nodes
std::shared_ptr< mesh::MeshNodes > mesh_nodes
Optional primitive-to-node mapping for this FE space.
Definition FESpace.hpp:86
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:437
polyfem::solver::AdjointTools::integrate_objective
double integrate_objective(const varform::DifferentiableVarForm &varform, 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::AdjointTools::dJ_shape_homogenization_adjoint_term
void dJ_shape_homogenization_adjoint_term(const varform::DifferentiableVarForm &varform, const DiffCache &diff_cache, const Eigen::MatrixXd &sol, const Eigen::MatrixXd &adjoint, Eigen::VectorXd &one_form)
Definition AdjointTools.cpp:636
polyfem::solver::Diff
DScalar1< double, Eigen::Matrix< double, Eigen::Dynamic, 1 > > Diff
Definition SurfaceTractionForms.cpp:251
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
polyfem::varform::VarFormBoundaryState::total_local_boundary
std::vector< mesh::LocalBoundary > total_local_boundary
Definition FESpace.hpp:154