PolyFEM
Loading...
Searching...
No Matches
AdjointTools.cpp
Go to the documentation of this file.
1#include "AdjointTools.hpp"
2
3#include <polyfem/utils/MaybeParallelFor.hpp>
4#include <polyfem/utils/Timer.hpp>
5#include <polyfem/io/Evaluator.hpp>
6#include <polyfem/utils/AutodiffTypes.hpp>
7#include <polyfem/State.hpp>
8
9#include <polyfem/utils/IntegrableFunctional.hpp>
10#include <polyfem/utils/BoundarySampler.hpp>
11
12#include <polyfem/solver/forms/ElasticForm.hpp>
13#include <polyfem/solver/forms/ContactForm.hpp>
14#include <polyfem/solver/forms/FrictionForm.hpp>
15#include <polyfem/solver/forms/BodyForm.hpp>
16#include <polyfem/solver/forms/PressureForm.hpp>
17#include <polyfem/solver/forms/InertiaForm.hpp>
18
19#include <polyfem/solver/NLProblem.hpp>
20
21#include <polyfem/time_integrator/BDF.hpp>
22
23/*
24Reminders:
25
26 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.
27
28*/
29
30using namespace polyfem::utils;
31
32namespace
33{
34
35 int get_bdf_order(const polyfem::State &state)
36 {
37 if (state.args["time"]["integrator"].is_string())
38 return 1;
39 if (state.args["time"]["integrator"]["type"] == "ImplicitEuler")
40 return 1;
41 if (state.args["time"]["integrator"]["type"] == "BDF")
42 return state.args["time"]["integrator"]["steps"].get<int>();
43
44 polyfem::log_and_throw_adjoint_error("Integrator type not supported for differentiability.");
45 return -1;
46 }
47} // namespace
48
49namespace polyfem::solver
50{
51 namespace
52 {
53 double dot(const Eigen::MatrixXd &A, const Eigen::MatrixXd &B) { return (A.array() * B.array()).sum(); }
54
55 class LocalThreadScalarStorage
56 {
57 public:
58 double val;
59 assembler::ElementAssemblyValues vals;
60 QuadratureVector da;
61
62 LocalThreadScalarStorage()
63 {
64 val = 0;
65 }
66 };
67
68 class LocalThreadVecStorage
69 {
70 public:
71 Eigen::MatrixXd vec;
72 assembler::ElementAssemblyValues vals;
73 QuadratureVector da;
74
75 LocalThreadVecStorage(const int size)
76 {
77 vec.resize(size, 1);
78 vec.setZero();
79 }
80 };
81
82 typedef DScalar1<double, Eigen::Matrix<double, Eigen::Dynamic, 1>> Diff;
83
84 template <typename T>
85 T triangle_area(const Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> &V)
86 {
87 Eigen::Matrix<T, Eigen::Dynamic, 1> l1 = V.row(1) - V.row(0);
88 Eigen::Matrix<T, Eigen::Dynamic, 1> l2 = V.row(2) - V.row(0);
89 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));
90 return area;
91 }
92
93 Eigen::MatrixXd triangle_area_grad(const Eigen::MatrixXd &F)
94 {
95 DiffScalarBase::setVariableCount(F.size());
96 Eigen::Matrix<Diff, Eigen::Dynamic, Eigen::Dynamic> full_diff(F.rows(), F.cols());
97 for (int i = 0; i < F.rows(); i++)
98 for (int j = 0; j < F.cols(); j++)
99 full_diff(i, j) = Diff(i + j * F.rows(), F(i, j));
100 auto reduced_diff = triangle_area(full_diff);
101
102 Eigen::MatrixXd grad(F.rows(), F.cols());
103 for (int i = 0; i < F.rows(); ++i)
104 for (int j = 0; j < F.cols(); ++j)
105 grad(i, j) = reduced_diff.getGradient()(i + j * F.rows());
106
107 return grad;
108 }
109
110 template <typename T>
111 T line_length(const Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> &V)
112 {
113 Eigen::Matrix<T, Eigen::Dynamic, 1> L = V.row(1) - V.row(0);
114 T area = L.norm();
115 return area;
116 }
117
118 Eigen::MatrixXd line_length_grad(const Eigen::MatrixXd &F)
119 {
120 DiffScalarBase::setVariableCount(F.size());
121 Eigen::Matrix<Diff, Eigen::Dynamic, Eigen::Dynamic> full_diff(F.rows(), F.cols());
122 for (int i = 0; i < F.rows(); i++)
123 for (int j = 0; j < F.cols(); j++)
124 full_diff(i, j) = Diff(i + j * F.rows(), F(i, j));
125 auto reduced_diff = line_length(full_diff);
126
127 Eigen::MatrixXd grad(F.rows(), F.cols());
128 for (int i = 0; i < F.rows(); ++i)
129 for (int j = 0; j < F.cols(); ++j)
130 grad(i, j) = reduced_diff.getGradient()(i + j * F.rows());
131
132 return grad;
133 }
134
135 template <typename T>
136 Eigen::Matrix<T, 2, 1> edge_normal(const Eigen::Matrix<T, 4, 1> &V)
137 {
138 Eigen::Matrix<T, 2, 1> v1 = V.segment(0, 2);
139 Eigen::Matrix<T, 2, 1> v2 = V.segment(2, 2);
140 Eigen::Matrix<T, 2, 1> normal = v1 - v2;
141 normal(0) *= -1;
142 normal = normal / normal.norm();
143 return normal;
144 }
145
146 template <typename T>
147 Eigen::Matrix<T, 3, 1> face_normal(const Eigen::Matrix<T, 9, 1> &V)
148 {
149 Eigen::Matrix<T, 3, 1> v1 = V.segment(0, 3);
150 Eigen::Matrix<T, 3, 1> v2 = V.segment(3, 3);
151 Eigen::Matrix<T, 3, 1> v3 = V.segment(6, 3);
152 Eigen::Matrix<T, 3, 1> normal = (v2 - v1).cross(v3 - v1);
153 normal = normal / normal.norm();
154 return normal;
155 }
156
157 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)
158 {
159 Eigen::MatrixXd params = Eigen::MatrixXd::Zero(local_pts.rows(), 2);
160
161 auto search_lambda = lame_params.find("lambda");
162 auto search_mu = lame_params.find("mu");
163
164 if (search_lambda == lame_params.end() || search_mu == lame_params.end())
165 return params;
166
167 for (int p = 0; p < local_pts.rows(); p++)
168 {
169 params(p, 0) = search_lambda->second(local_pts.row(p), pts.row(p), t, e);
170 params(p, 1) = search_mu->second(local_pts.row(p), pts.row(p), t, e);
171 }
172
173 return params;
174 }
175 } // namespace
176
177 double AdjointTools::integrate_objective(
178 const State &state,
179 const IntegrableFunctional &j,
180 const Eigen::MatrixXd &solution,
181 const std::set<int> &interested_ids, // either body id or surface id
182 const SpatialIntegralType spatial_integral_type,
183 const int cur_step) // current time step
184 {
185 const auto &bases = state.bases;
186 const auto &gbases = state.geom_bases();
187
188 const int dim = state.mesh->dimension();
189 const int actual_dim = state.problem->is_scalar() ? 1 : dim;
190 const int n_elements = int(bases.size());
191 const double t0 = state.problem->is_time_dependent() ? state.args["time"]["t0"].get<double>() : 0.0;
192 const double dt = state.problem->is_time_dependent() ? state.args["time"]["dt"].get<double>() : 0.0;
193
194 double integral = 0;
195 if (spatial_integral_type == SpatialIntegralType::Volume)
196 {
197 auto storage = utils::create_thread_storage(LocalThreadScalarStorage());
198 utils::maybe_parallel_for(n_elements, [&](int start, int end, int thread_id) {
199 LocalThreadScalarStorage &local_storage = utils::get_local_thread_storage(storage, thread_id);
200
201 IntegrableFunctional::ParameterType params;
202 params.t = dt * cur_step + t0;
203 params.step = cur_step;
204
205 Eigen::MatrixXd u, grad_u;
206 Eigen::MatrixXd result;
207
208 for (int e = start; e < end; ++e)
209 {
210 if (interested_ids.size() != 0 && interested_ids.find(state.mesh->get_body_id(e)) == interested_ids.end())
211 continue;
212
213 assembler::ElementAssemblyValues &vals = local_storage.vals;
214 state.ass_vals_cache.compute(e, state.mesh->is_volume(), bases[e], gbases[e], vals);
215 io::Evaluator::interpolate_at_local_vals(e, dim, actual_dim, vals, solution, u, grad_u);
216
217 const quadrature::Quadrature &quadrature = vals.quadrature;
218 local_storage.da = vals.det.array() * quadrature.weights.array();
219
220 const Eigen::MatrixXd lame_params = extract_lame_params(state.assembler->parameters(), e, params.t, quadrature.points, vals.val);
221
222 params.elem = e;
223 params.body_id = state.mesh->get_body_id(e);
224 j.evaluate(lame_params, quadrature.points, vals.val, u, grad_u, Eigen::MatrixXd::Zero(0, 0) /*Not used*/, vals, params, result);
225
226 local_storage.val += dot(result, local_storage.da);
227 }
228 });
229 for (const LocalThreadScalarStorage &local_storage : storage)
230 integral += local_storage.val;
231 }
232 else if (spatial_integral_type == SpatialIntegralType::Surface)
233 {
234 auto storage = utils::create_thread_storage(LocalThreadScalarStorage());
235 utils::maybe_parallel_for(state.total_local_boundary.size(), [&](int start, int end, int thread_id) {
236 LocalThreadScalarStorage &local_storage = utils::get_local_thread_storage(storage, thread_id);
237
238 Eigen::MatrixXd uv;
239 Eigen::MatrixXd points, normal;
240 Eigen::VectorXd weights;
241
242 Eigen::MatrixXd u, grad_u;
243 Eigen::MatrixXd result;
244 IntegrableFunctional::ParameterType params;
245 params.t = dt * cur_step + t0;
246 params.step = cur_step;
247
248 for (int lb_id = start; lb_id < end; ++lb_id)
249 {
250 const auto &lb = state.total_local_boundary[lb_id];
251 const int e = lb.element_id();
252
253 for (int i = 0; i < lb.size(); i++)
254 {
255 const int global_primitive_id = lb.global_primitive_id(i);
256 if (interested_ids.size() != 0 && interested_ids.find(state.mesh->get_boundary_id(global_primitive_id)) == interested_ids.end())
257 continue;
258
259 utils::BoundarySampler::boundary_quadrature(lb, state.n_boundary_samples(), *state.mesh, i, false, uv, points, normal, weights);
260
261 assembler::ElementAssemblyValues &vals = local_storage.vals;
262 vals.compute(e, state.mesh->is_volume(), points, bases[e], gbases[e]);
263 io::Evaluator::interpolate_at_local_vals(e, dim, actual_dim, vals, solution, u, grad_u);
264
265 const Eigen::MatrixXd lame_params = extract_lame_params(state.assembler->parameters(), e, params.t, points, vals.val);
266
267 params.elem = e;
268 params.body_id = state.mesh->get_body_id(e);
269 params.boundary_id = state.mesh->get_boundary_id(global_primitive_id);
270 j.evaluate(lame_params, points, vals.val, u, grad_u, normal, vals, params, result);
271
272 local_storage.val += dot(result, weights);
273 }
274 }
275 });
276 for (const LocalThreadScalarStorage &local_storage : storage)
277 integral += local_storage.val;
278 }
279 else if (spatial_integral_type == SpatialIntegralType::VertexSum)
280 {
281 std::vector<bool> traversed(state.n_bases, false);
282 IntegrableFunctional::ParameterType params;
283 params.t = dt * cur_step + t0;
284 params.step = cur_step;
285 for (int e = 0; e < bases.size(); e++)
286 {
287 const auto &bs = bases[e];
288 for (int i = 0; i < bs.bases.size(); i++)
289 {
290 const auto &b = bs.bases[i];
291 assert(b.global().size() == 1);
292 const auto &g = b.global()[0];
293 if (traversed[g.index])
294 continue;
295
296 const Eigen::MatrixXd lame_params = extract_lame_params(state.assembler->parameters(), e, params.t, Eigen::MatrixXd::Zero(1, dim) /*Not used*/, g.node);
297
298 params.node = g.index;
299 params.elem = e;
300 params.body_id = state.mesh->get_body_id(e);
301 Eigen::MatrixXd val;
302 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);
303 integral += val(0);
304 traversed[g.index] = true;
305 }
306 }
307 }
308
309 return integral;
310 }
311
312 void AdjointTools::compute_shape_derivative_functional_term(
313 const State &state,
314 const Eigen::MatrixXd &solution,
315 const IntegrableFunctional &j,
316 const std::set<int> &interested_ids, // either body id or surface id
317 const SpatialIntegralType spatial_integral_type,
318 Eigen::VectorXd &term,
319 const int cur_time_step)
320 {
321 const auto &gbases = state.geom_bases();
322 const auto &bases = state.bases;
323 const int dim = state.mesh->dimension();
324 const int actual_dim = state.problem->is_scalar() ? 1 : dim;
325 const double t0 = state.problem->is_time_dependent() ? state.args["time"]["t0"].get<double>() : 0.0;
326 const double dt = state.problem->is_time_dependent() ? state.args["time"]["dt"].get<double>() : 0.0;
327
328 const int n_elements = int(bases.size());
329 term.setZero(state.n_geom_bases * dim, 1);
330
331 auto storage = utils::create_thread_storage(LocalThreadVecStorage(term.size()));
332
333 if (spatial_integral_type == SpatialIntegralType::Volume)
334 {
335 utils::maybe_parallel_for(n_elements, [&](int start, int end, int thread_id) {
336 LocalThreadVecStorage &local_storage = utils::get_local_thread_storage(storage, thread_id);
337
338 Eigen::MatrixXd u, grad_u, j_val, dj_dgradu, dj_dx;
339
340 IntegrableFunctional::ParameterType params;
341 params.t = cur_time_step * dt + t0;
342 params.step = cur_time_step;
343
344 for (int e = start; e < end; ++e)
345 {
346 if (interested_ids.size() != 0 && interested_ids.find(state.mesh->get_body_id(e)) == interested_ids.end())
347 continue;
348
349 assembler::ElementAssemblyValues &vals = local_storage.vals;
350 state.ass_vals_cache.compute(e, state.mesh->is_volume(), bases[e], gbases[e], vals);
351 io::Evaluator::interpolate_at_local_vals(e, dim, actual_dim, vals, solution, u, grad_u);
352
353 assembler::ElementAssemblyValues gvals;
354 gvals.compute(e, state.mesh->is_volume(), vals.quadrature.points, gbases[e], gbases[e]);
355
356 const quadrature::Quadrature &quadrature = vals.quadrature;
357 local_storage.da = vals.det.array() * quadrature.weights.array();
358
359 const Eigen::MatrixXd lame_params = extract_lame_params(state.assembler->parameters(), e, params.t, quadrature.points, vals.val);
360
361 params.elem = e;
362 params.body_id = state.mesh->get_body_id(e);
363
364 j.evaluate(lame_params, quadrature.points, vals.val, u, grad_u, Eigen::MatrixXd::Zero(0, 0) /*Not used*/, vals, params, j_val);
365
366 if (j.depend_on_gradu())
367 j.dj_dgradu(lame_params, quadrature.points, vals.val, u, grad_u, Eigen::MatrixXd::Zero(0, 0) /*Not used*/, vals, params, dj_dgradu);
368
369 if (j.depend_on_x())
370 j.dj_dx(lame_params, quadrature.points, vals.val, u, grad_u, Eigen::MatrixXd::Zero(0, 0) /*Not used*/, vals, params, dj_dx);
371
372 Eigen::MatrixXd tau_q, grad_u_q;
373 for (auto &v : gvals.basis_values)
374 {
375 for (int q = 0; q < local_storage.da.size(); ++q)
376 {
377 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();
378
379 if (j.depend_on_x())
380 local_storage.vec.block(v.global[0].index * dim, 0, dim, 1) += (v.val(q) * local_storage.da(q)) * dj_dx.row(q).transpose();
381
382 if (j.depend_on_gradu())
383 {
384 if (dim == actual_dim) // Elasticity PDE
385 {
386 vector2matrix(dj_dgradu.row(q), tau_q);
387 vector2matrix(grad_u.row(q), grad_u_q);
388 }
389 else // Laplacian PDE
390 {
391 tau_q = dj_dgradu.row(q);
392 grad_u_q = grad_u.row(q);
393 }
394 for (int d = 0; d < dim; d++)
395 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);
396 }
397 }
398 }
399 }
400 });
401 }
402 else if (spatial_integral_type == SpatialIntegralType::Surface)
403 {
404 utils::maybe_parallel_for(state.total_local_boundary.size(), [&](int start, int end, int thread_id) {
405 LocalThreadVecStorage &local_storage = utils::get_local_thread_storage(storage, thread_id);
406
407 Eigen::MatrixXd uv, points, normal;
408 Eigen::VectorXd &weights = local_storage.da;
409
410 Eigen::MatrixXd u, grad_u, x, grad_x, j_val, dj_dgradu, dj_dgradx, dj_dx;
411
412 IntegrableFunctional::ParameterType params;
413 params.t = cur_time_step * dt + t0;
414 params.step = cur_time_step;
415
416 for (int lb_id = start; lb_id < end; ++lb_id)
417 {
418 const auto &lb = state.total_local_boundary[lb_id];
419 const int e = lb.element_id();
420
421 for (int i = 0; i < lb.size(); i++)
422 {
423 const int global_primitive_id = lb.global_primitive_id(i);
424 if (interested_ids.size() != 0 && interested_ids.find(state.mesh->get_boundary_id(global_primitive_id)) == interested_ids.end())
425 continue;
426
427 utils::BoundarySampler::boundary_quadrature(lb, state.n_boundary_samples(), *state.mesh, i, false, uv, points, normal, weights);
428
429 assembler::ElementAssemblyValues &vals = local_storage.vals;
430 io::Evaluator::interpolate_at_local_vals(*state.mesh, state.problem->is_scalar(), bases, gbases, e, points, solution, u, grad_u);
431 // io::Evaluator::interpolate_at_local_vals(*state.mesh, state.problem->is_scalar(), gbases, gbases, e, points, global_positions, x, grad_x);
432
433 vals.compute(e, state.mesh->is_volume(), points, gbases[e], gbases[e]);
434
435 // normal = normal * vals.jac_it[0]; // assuming linear geometry
436
437 const int n_loc_bases_ = int(vals.basis_values.size());
438
439 const Eigen::MatrixXd lame_params = extract_lame_params(state.assembler->parameters(), e, params.t, points, vals.val);
440
441 params.elem = e;
442 params.body_id = state.mesh->get_body_id(e);
443 params.boundary_id = state.mesh->get_boundary_id(global_primitive_id);
444
445 j.evaluate(lame_params, points, vals.val, u, grad_u, normal, vals, params, j_val);
446 j_val = j_val.array().colwise() * weights.array();
447
448 if (j.depend_on_gradu())
449 {
450 j.dj_dgradu(lame_params, points, vals.val, u, grad_u, normal, vals, params, dj_dgradu);
451 dj_dgradu = dj_dgradu.array().colwise() * weights.array();
452 }
453
454 if (j.depend_on_gradx())
455 {
456 j.dj_dgradx(lame_params, points, vals.val, u, grad_u, normal, vals, params, dj_dgradx);
457 dj_dgradx = dj_dgradx.array().colwise() * weights.array();
458 }
459
460 if (j.depend_on_x())
461 {
462 j.dj_dx(lame_params, points, vals.val, u, grad_u, normal, vals, params, dj_dx);
463 dj_dx = dj_dx.array().colwise() * weights.array();
464 }
465
466 const auto nodes = gbases[e].local_nodes_for_primitive(lb.global_primitive_id(i), *state.mesh);
467
468 if (nodes.size() != dim)
469 log_and_throw_adjoint_error("Only linear geometry is supported in differentiable surface integral functional!");
470
471 Eigen::MatrixXd velocity_div_mat;
472 if (state.mesh->is_volume())
473 {
474 Eigen::Matrix3d V;
475 for (int d = 0; d < 3; d++)
476 V.row(d) = gbases[e].bases[nodes(d)].global()[0].node;
477 velocity_div_mat = face_velocity_divergence(V);
478 }
479 else
480 {
481 Eigen::Matrix2d V;
482 for (int d = 0; d < 2; d++)
483 V.row(d) = gbases[e].bases[nodes(d)].global()[0].node;
484 velocity_div_mat = edge_velocity_divergence(V);
485 }
486
487 Eigen::MatrixXd grad_u_q, tau_q, grad_x_q;
488 for (long n = 0; n < nodes.size(); ++n)
489 {
490 const assembler::AssemblyValues &v = vals.basis_values[nodes(n)];
491
492 local_storage.vec.block(v.global[0].index * dim, 0, dim, 1) += j_val.sum() * velocity_div_mat.row(n).transpose();
493 }
494
495 for (long n = 0; n < n_loc_bases_; ++n)
496 {
497 const assembler::AssemblyValues &v = vals.basis_values[n];
498
499 if (j.depend_on_x())
500 local_storage.vec.block(v.global[0].index * dim, 0, dim, 1) += dj_dx.transpose() * v.val;
501
502 // integrate j * div(gbases) over the whole boundary
503 if (j.depend_on_gradu())
504 {
505 for (int q = 0; q < weights.size(); ++q)
506 {
507 if (dim == actual_dim) // Elasticity PDE
508 {
509 vector2matrix(grad_u.row(q), grad_u_q);
510 vector2matrix(dj_dgradu.row(q), tau_q);
511 }
512 else // Laplacian PDE
513 {
514 grad_u_q = grad_u.row(q);
515 tau_q = dj_dgradu.row(q);
516 }
517
518 for (int d = 0; d < dim; d++)
519 local_storage.vec(v.global[0].index * dim + d) += -dot(tau_q, grad_u_q.col(d) * v.grad_t_m.row(q));
520 }
521 }
522
523 if (j.depend_on_gradx())
524 {
525 for (int d = 0; d < dim; d++)
526 {
527 for (int q = 0; q < weights.size(); ++q)
528 local_storage.vec(v.global[0].index * dim + d) += dot(dj_dgradx.block(q, d * dim, 1, dim), v.grad.row(q));
529 }
530 }
531 }
532 }
533 }
534 });
535 }
536 else if (spatial_integral_type == SpatialIntegralType::VertexSum)
537 {
538 log_and_throw_adjoint_error("Shape derivative of vertex sum type functional is not implemented!");
539 }
540 for (const LocalThreadVecStorage &local_storage : storage)
541 term += local_storage.vec;
542
543 term = utils::flatten(utils::unflatten(term, dim)(state.primitive_to_node(), Eigen::all));
544 }
545
546 void AdjointTools::dJ_shape_static_adjoint_term(
547 const State &state,
548 const Eigen::MatrixXd &sol,
549 const Eigen::MatrixXd &adjoint,
550 Eigen::VectorXd &one_form)
551 {
552 Eigen::VectorXd elasticity_term, rhs_term, pressure_term, contact_term;
553
554 one_form.setZero(state.n_geom_bases * state.mesh->dimension());
555
556 // if (j.depend_on_u() || j.depend_on_gradu())
557 {
558 state.solve_data.elastic_form->force_shape_derivative(0, state.n_geom_bases, sol, sol, adjoint, elasticity_term);
559 if (state.solve_data.body_form)
560 state.solve_data.body_form->force_shape_derivative(state.n_geom_bases, 0, sol, adjoint, rhs_term);
561 else
562 rhs_term.setZero(one_form.size());
563
564 if (state.solve_data.pressure_form)
565 {
566 state.solve_data.pressure_form->force_shape_derivative(state.n_geom_bases, 0, sol, adjoint, pressure_term);
567 pressure_term = state.gbasis_nodes_to_basis_nodes * pressure_term;
568 }
569 else
570 pressure_term.setZero(one_form.size());
571
572 if (state.is_contact_enabled())
573 {
574 state.solve_data.contact_form->force_shape_derivative(state.diff_cached.collision_set(0), sol, adjoint, contact_term);
575 contact_term = state.gbasis_nodes_to_basis_nodes * contact_term;
576 }
577 else
578 contact_term.setZero(elasticity_term.size());
579 one_form -= elasticity_term + rhs_term + pressure_term + contact_term;
580 }
581
582 one_form = utils::flatten(utils::unflatten(one_form, state.mesh->dimension())(state.primitive_to_node(), Eigen::all));
583 }
584
585 void AdjointTools::dJ_shape_transient_adjoint_term(
586 const State &state,
587 const Eigen::MatrixXd &adjoint_nu,
588 const Eigen::MatrixXd &adjoint_p,
589 Eigen::VectorXd &one_form)
590 {
591 const double t0 = state.args["time"]["t0"];
592 const double dt = state.args["time"]["dt"];
593 const int time_steps = state.args["time"]["time_steps"];
594 const int bdf_order = get_bdf_order(state);
595
596 Eigen::VectorXd elasticity_term, rhs_term, pressure_term, damping_term, mass_term, contact_term, friction_term;
597 one_form.setZero(state.n_geom_bases * state.mesh->dimension());
598
599 Eigen::VectorXd cur_p, cur_nu;
600 for (int i = time_steps; i > 0; --i)
601 {
602 const int real_order = std::min(bdf_order, i);
603 double beta = time_integrator::BDF::betas(real_order - 1);
604 double beta_dt = beta * dt;
605 const double t = i * dt + t0;
606
607 Eigen::MatrixXd velocity = state.diff_cached.v(i);
608
609 cur_p = adjoint_p.col(i);
610 cur_nu = adjoint_nu.col(i);
611 cur_p(state.boundary_nodes).setZero();
612 cur_nu(state.boundary_nodes).setZero();
613
614 {
615 state.solve_data.inertia_form->force_shape_derivative(state.mesh->is_volume(), state.n_geom_bases, t, state.bases, state.geom_bases(), *(state.mass_matrix_assembler), state.mass_ass_vals_cache, velocity, cur_nu, mass_term);
616 state.solve_data.elastic_form->force_shape_derivative(t, state.n_geom_bases, state.diff_cached.u(i), state.diff_cached.u(i), cur_p, elasticity_term);
617 state.solve_data.body_form->force_shape_derivative(state.n_geom_bases, t, state.diff_cached.u(i - 1), cur_p, rhs_term);
618 state.solve_data.pressure_form->force_shape_derivative(state.n_geom_bases, t, state.diff_cached.u(i), cur_p, pressure_term);
619 pressure_term = state.gbasis_nodes_to_basis_nodes * pressure_term;
620
621 if (state.solve_data.damping_form)
622 state.solve_data.damping_form->force_shape_derivative(t, state.n_geom_bases, state.diff_cached.u(i), state.diff_cached.u(i - 1), cur_p, damping_term);
623 else
624 damping_term.setZero(mass_term.size());
625
626 if (state.is_contact_enabled())
627 {
628 state.solve_data.contact_form->force_shape_derivative(state.diff_cached.collision_set(i), state.diff_cached.u(i), cur_p, contact_term);
629 contact_term = state.gbasis_nodes_to_basis_nodes * contact_term;
630 // contact_term /= beta_dt * beta_dt;
631 }
632 else
633 contact_term.setZero(mass_term.size());
634
635 if (state.solve_data.friction_form)
636 {
637 state.solve_data.friction_form->force_shape_derivative(state.diff_cached.u(i - 1), state.diff_cached.u(i), cur_p, state.diff_cached.friction_collision_set(i), friction_term);
638 friction_term = state.gbasis_nodes_to_basis_nodes * (friction_term / beta);
639 // friction_term /= beta_dt * beta_dt;
640 }
641 else
642 friction_term.setZero(mass_term.size());
643 }
644
645 one_form += beta_dt * (elasticity_term + rhs_term + pressure_term + damping_term + contact_term + friction_term + mass_term);
646 }
647
648 // time step 0
649 Eigen::VectorXd sum_alpha_p;
650 {
651 sum_alpha_p.setZero(adjoint_p.rows());
652 int num = std::min(bdf_order, time_steps);
653 for (int j = 0; j < num; ++j)
654 {
655 int order = std::min(bdf_order - 1, j);
656 sum_alpha_p -= time_integrator::BDF::alphas(order)[j] * adjoint_p.col(j + 1);
657 }
658 }
659 sum_alpha_p(state.boundary_nodes).setZero();
660 state.solve_data.inertia_form->force_shape_derivative(state.mesh->is_volume(), state.n_geom_bases, t0, state.bases, state.geom_bases(), *(state.mass_matrix_assembler), state.mass_ass_vals_cache, state.diff_cached.v(0), sum_alpha_p, mass_term);
661
662 one_form += mass_term;
663
664 one_form = utils::flatten(utils::unflatten(one_form, state.mesh->dimension())(state.primitive_to_node(), Eigen::all));
665 }
666
667 void AdjointTools::dJ_material_static_adjoint_term(
668 const State &state,
669 const Eigen::MatrixXd &sol,
670 const Eigen::MatrixXd &adjoint,
671 Eigen::VectorXd &one_form)
672 {
673 state.solve_data.elastic_form->force_material_derivative(0, sol, sol, adjoint, one_form);
674 }
675
676 void AdjointTools::dJ_material_transient_adjoint_term(
677 const State &state,
678 const Eigen::MatrixXd &adjoint_nu,
679 const Eigen::MatrixXd &adjoint_p,
680 Eigen::VectorXd &one_form)
681 {
682 const double t0 = state.args["time"]["t0"];
683 const double dt = state.args["time"]["dt"];
684 const int time_steps = state.args["time"]["time_steps"];
685 const int bdf_order = get_bdf_order(state);
686
687 one_form.setZero(state.bases.size() * 2);
688
689 auto storage = utils::create_thread_storage(LocalThreadVecStorage(one_form.size()));
690
691 utils::maybe_parallel_for(time_steps, [&](int start, int end, int thread_id) {
692 LocalThreadVecStorage &local_storage = utils::get_local_thread_storage(storage, thread_id);
693 Eigen::VectorXd elasticity_term;
694 for (int i_aux = start; i_aux < end; ++i_aux)
695 {
696 const int i = time_steps - i_aux;
697 const int real_order = std::min(bdf_order, i);
698 double beta_dt = time_integrator::BDF::betas(real_order - 1) * dt;
699
700 Eigen::VectorXd cur_p = adjoint_p.col(i);
701 cur_p(state.boundary_nodes).setZero();
702
703 state.solve_data.elastic_form->force_material_derivative(t0 + dt * i, state.diff_cached.u(i), state.diff_cached.u(i - 1), -cur_p, elasticity_term);
704 local_storage.vec += beta_dt * elasticity_term;
705 }
706 });
707
708 for (const LocalThreadVecStorage &local_storage : storage)
709 one_form += local_storage.vec;
710 }
711
712 void AdjointTools::dJ_friction_transient_adjoint_term(
713 const State &state,
714 const Eigen::MatrixXd &adjoint_nu,
715 const Eigen::MatrixXd &adjoint_p,
716 Eigen::VectorXd &one_form)
717 {
718 const double dt = state.args["time"]["dt"];
719 const double mu = state.solve_data.friction_form->mu();
720 const int time_steps = state.args["time"]["time_steps"];
721 const int dim = state.mesh->dimension();
722 const int bdf_order = get_bdf_order(state);
723
724 one_form.setZero(1);
725
726 auto storage = utils::create_thread_storage(LocalThreadScalarStorage());
727
728 utils::maybe_parallel_for(time_steps, [&](int start, int end, int thread_id) {
729 LocalThreadScalarStorage &local_storage = utils::get_local_thread_storage(storage, thread_id);
730 for (int t_aux = start; t_aux < end; ++t_aux)
731 {
732 const int t = time_steps - t_aux;
733 const int real_order = std::min(bdf_order, t);
734 double beta = time_integrator::BDF::betas(real_order - 1);
735
736 const Eigen::MatrixXd surface_solution_prev = state.collision_mesh.vertices(utils::unflatten(state.diff_cached.u(t - 1), dim));
737 const Eigen::MatrixXd surface_solution = state.collision_mesh.vertices(utils::unflatten(state.diff_cached.u(t), dim));
738
739 // TODO: use the time integration to compute the velocity
740 const Eigen::MatrixXd surface_velocities = (surface_solution - surface_solution_prev) / dt;
741
742 Eigen::MatrixXd force = state.collision_mesh.to_full_dof(
743 -state.solve_data.friction_form->friction_potential().force(
744 state.diff_cached.friction_collision_set(t),
745 state.collision_mesh,
746 state.collision_mesh.rest_positions(),
747 /*lagged_displacements=*/surface_solution_prev,
748 surface_velocities,
749 state.solve_data.contact_form->barrier_potential(),
750 state.solve_data.contact_form->barrier_stiffness(),
751 state.solve_data.friction_form->epsv()));
752
753 Eigen::VectorXd cur_p = adjoint_p.col(t);
754 cur_p(state.boundary_nodes).setZero();
755
756 local_storage.val += dot(cur_p, force) / (beta * mu * dt);
757 }
758 });
759
760 for (const LocalThreadScalarStorage &local_storage : storage)
761 one_form(0) += local_storage.val;
762 }
763
764 void AdjointTools::dJ_damping_transient_adjoint_term(
765 const State &state,
766 const Eigen::MatrixXd &adjoint_nu,
767 const Eigen::MatrixXd &adjoint_p,
768 Eigen::VectorXd &one_form)
769 {
770 const double t0 = state.args["time"]["t0"];
771 const double dt = state.args["time"]["dt"];
772 const int time_steps = state.args["time"]["time_steps"];
773 const int bdf_order = get_bdf_order(state);
774
775 one_form.setZero(2);
776
777 auto storage = utils::create_thread_storage(LocalThreadVecStorage(one_form.size()));
778
779 utils::maybe_parallel_for(time_steps, [&](int start, int end, int thread_id) {
780 LocalThreadVecStorage &local_storage = utils::get_local_thread_storage(storage, thread_id);
781 Eigen::VectorXd damping_term;
782 for (int t_aux = start; t_aux < end; ++t_aux)
783 {
784 const int t = time_steps - t_aux;
785 const int real_order = std::min(bdf_order, t);
786 const double beta = time_integrator::BDF::betas(real_order - 1);
787
788 Eigen::VectorXd cur_p = adjoint_p.col(t);
789 cur_p(state.boundary_nodes).setZero();
790
791 state.solve_data.damping_form->force_material_derivative(t * dt + t0, state.diff_cached.u(t), state.diff_cached.u(t - 1), -cur_p, damping_term);
792 local_storage.vec += (beta * dt) * damping_term;
793 }
794 });
795
796 for (const LocalThreadVecStorage &local_storage : storage)
797 one_form += local_storage.vec;
798 }
799
800 void AdjointTools::dJ_initial_condition_adjoint_term(
801 const State &state,
802 const Eigen::MatrixXd &adjoint_nu,
803 const Eigen::MatrixXd &adjoint_p,
804 Eigen::VectorXd &one_form)
805 {
806 const int ndof = state.ndof();
807 one_form.setZero(ndof * 2); // half for initial solution, half for initial velocity
808
809 // \partial_q \hat{J}^0 - p_0^T \partial_q g^v - \mu_0^T \partial_q g^u
810 one_form.segment(0, ndof) = -adjoint_nu.col(0); // adjoint_nu[0] actually stores adjoint_mu[0]
811 one_form.segment(ndof, ndof) = -adjoint_p.col(0);
812
813 for (int b : state.boundary_nodes)
814 {
815 one_form(b) = 0;
816 one_form(ndof + b) = 0;
817 }
818 }
819
820 void AdjointTools::dJ_dirichlet_transient_adjoint_term(
821 const State &state,
822 const Eigen::MatrixXd &adjoint_nu,
823 const Eigen::MatrixXd &adjoint_p,
824 Eigen::VectorXd &one_form)
825 {
826 const double dt = state.args["time"]["dt"];
827 const int time_steps = state.args["time"]["time_steps"];
828 const int bdf_order = get_bdf_order(state);
829 const int n_dirichlet_dof = state.boundary_nodes.size();
830
831 one_form.setZero(time_steps * n_dirichlet_dof);
832 for (int i = 1; i <= time_steps; ++i)
833 {
834 const int real_order = std::min(bdf_order, i);
835 const double beta_dt = time_integrator::BDF::betas(real_order - 1) * dt;
836
837 one_form.segment((i - 1) * n_dirichlet_dof, n_dirichlet_dof) = -(1. / beta_dt) * adjoint_p(state.boundary_nodes, i);
838 }
839 }
840
841 void AdjointTools::dJ_du_step(
842 const State &state,
843 const IntegrableFunctional &j,
844 const Eigen::MatrixXd &solution,
845 const std::set<int> &interested_ids,
846 const SpatialIntegralType spatial_integral_type,
847 const int cur_step,
848 Eigen::VectorXd &term)
849 {
850 const auto &bases = state.bases;
851 const auto &gbases = state.geom_bases();
852
853 const int dim = state.mesh->dimension();
854 const int actual_dim = state.problem->is_scalar() ? 1 : dim;
855 const int n_elements = int(bases.size());
856 const double t0 = state.problem->is_time_dependent() ? state.args["time"]["t0"].get<double>() : 0.0;
857 const double dt = state.problem->is_time_dependent() ? state.args["time"]["dt"].get<double>() : 0.0;
858
859 term = Eigen::MatrixXd::Zero(state.n_bases * actual_dim, 1);
860
861 if (!j.depend_on_u() && !j.depend_on_gradu() && !j.depend_on_gradu_local())
862 return;
863
864 if (spatial_integral_type == SpatialIntegralType::Volume)
865 {
866 auto storage = utils::create_thread_storage(LocalThreadVecStorage(term.size()));
867 utils::maybe_parallel_for(n_elements, [&](int start, int end, int thread_id) {
868 LocalThreadVecStorage &local_storage = utils::get_local_thread_storage(storage, thread_id);
869
870 Eigen::MatrixXd u, grad_u;
871 Eigen::MatrixXd lambda, mu;
872 Eigen::MatrixXd dj_du, dj_dgradu, dj_dgradx;
873
874 IntegrableFunctional::ParameterType params;
875 params.t = dt * cur_step + t0;
876 params.step = cur_step;
877
878 for (int e = start; e < end; ++e)
879 {
880 if (interested_ids.size() != 0 && interested_ids.find(state.mesh->get_body_id(e)) == interested_ids.end())
881 continue;
882
883 assembler::ElementAssemblyValues &vals = local_storage.vals;
884 state.ass_vals_cache.compute(e, state.mesh->is_volume(), bases[e], gbases[e], vals);
885
886 const quadrature::Quadrature &quadrature = vals.quadrature;
887 local_storage.da = vals.det.array() * quadrature.weights.array();
888
889 const Eigen::MatrixXd lame_params = extract_lame_params(state.assembler->parameters(), e, params.t, quadrature.points, vals.val);
890
891 const int n_loc_bases_ = int(vals.basis_values.size());
892
893 io::Evaluator::interpolate_at_local_vals(e, dim, actual_dim, vals, solution, u, grad_u);
894
895 params.elem = e;
896 params.body_id = state.mesh->get_body_id(e);
897
898 dj_dgradu.resize(0, 0);
899 if (j.depend_on_gradu())
900 {
901 j.dj_dgradu(lame_params, quadrature.points, vals.val, u, grad_u, Eigen::MatrixXd::Zero(0, 0) /*Not used*/, vals, params, dj_dgradu);
902 for (int q = 0; q < dj_dgradu.rows(); q++)
903 dj_dgradu.row(q) *= local_storage.da(q);
904 }
905
906 dj_du.resize(0, 0);
907 if (j.depend_on_u())
908 {
909 j.dj_du(lame_params, quadrature.points, vals.val, u, grad_u, Eigen::MatrixXd::Zero(0, 0) /*Not used*/, vals, params, dj_du);
910 for (int q = 0; q < dj_du.rows(); q++)
911 dj_du.row(q) *= local_storage.da(q);
912 }
913
914 for (int i = 0; i < n_loc_bases_; ++i)
915 {
916 const assembler::AssemblyValues &v = vals.basis_values[i];
917 assert(v.global.size() == 1);
918 for (int d = 0; d < actual_dim; d++)
919 {
920 double val = 0;
921
922 // j = j(..., grad u)
923 if (j.depend_on_gradu())
924 {
925 for (int q = 0; q < local_storage.da.size(); ++q)
926 val += dot(dj_dgradu.block(q, d * dim, 1, dim), v.grad_t_m.row(q));
927 }
928
929 // j = j(..., u)
930 if (j.depend_on_u())
931 {
932 for (int q = 0; q < local_storage.da.size(); ++q)
933 val += dj_du(q, d) * v.val(q);
934 }
935 local_storage.vec(v.global[0].index * actual_dim + d) += val;
936 }
937 }
938 }
939 });
940 for (const LocalThreadVecStorage &local_storage : storage)
941 term += local_storage.vec;
942 }
943 else if (spatial_integral_type == SpatialIntegralType::Surface)
944 {
945 auto storage = utils::create_thread_storage(LocalThreadVecStorage(term.size()));
946 utils::maybe_parallel_for(state.total_local_boundary.size(), [&](int start, int end, int thread_id) {
947 LocalThreadVecStorage &local_storage = utils::get_local_thread_storage(storage, thread_id);
948
949 Eigen::MatrixXd uv, samples, gtmp;
950 Eigen::MatrixXd points, normal;
951 Eigen::VectorXd weights;
952
953 Eigen::MatrixXd u, grad_u;
954 Eigen::MatrixXd lambda, mu;
955 Eigen::MatrixXd dj_du, dj_dgradu, dj_dgradu_local;
956
957 IntegrableFunctional::ParameterType params;
958 params.t = dt * cur_step + t0;
959 params.step = cur_step;
960
961 for (int lb_id = start; lb_id < end; ++lb_id)
962 {
963 const auto &lb = state.total_local_boundary[lb_id];
964 const int e = lb.element_id();
965
966 for (int i = 0; i < lb.size(); i++)
967 {
968 const int global_primitive_id = lb.global_primitive_id(i);
969 if (interested_ids.size() != 0 && interested_ids.find(state.mesh->get_boundary_id(global_primitive_id)) == interested_ids.end())
970 continue;
971
972 utils::BoundarySampler::boundary_quadrature(lb, state.n_boundary_samples(), *state.mesh, i, false, uv, points, normal, weights);
973
974 assembler::ElementAssemblyValues &vals = local_storage.vals;
975 vals.compute(e, state.mesh->is_volume(), points, bases[e], gbases[e]);
976 io::Evaluator::interpolate_at_local_vals(e, dim, actual_dim, vals, solution, u, grad_u);
977
978 const Eigen::MatrixXd lame_params = extract_lame_params(state.assembler->parameters(), e, params.t, points, vals.val);
979
980 // normal = normal * vals.jac_it[0]; // assuming linear geometry
981
982 const int n_loc_bases_ = int(vals.basis_values.size());
983
984 params.elem = e;
985 params.body_id = state.mesh->get_body_id(e);
986 params.boundary_id = state.mesh->get_boundary_id(global_primitive_id);
987
988 dj_dgradu.resize(0, 0);
989 if (j.depend_on_gradu())
990 {
991 j.dj_dgradu(lame_params, points, vals.val, u, grad_u, normal, vals, params, dj_dgradu);
992 for (int q = 0; q < dj_dgradu.rows(); q++)
993 dj_dgradu.row(q) *= weights(q);
994 }
995
996 dj_dgradu_local.resize(0, 0);
997 if (j.depend_on_gradu_local())
998 {
999 j.dj_dgradu_local(lame_params, points, vals.val, u, grad_u, normal, vals, params, dj_dgradu_local);
1000 for (int q = 0; q < dj_dgradu_local.rows(); q++)
1001 dj_dgradu_local.row(q) *= weights(q);
1002 }
1003
1004 dj_du.resize(0, 0);
1005 if (j.depend_on_u())
1006 {
1007 j.dj_du(lame_params, points, vals.val, u, grad_u, normal, vals, params, dj_du);
1008 for (int q = 0; q < dj_du.rows(); q++)
1009 dj_du.row(q) *= weights(q);
1010 }
1011
1012 for (int l = 0; l < lb.size(); ++l)
1013 {
1014 const auto nodes = bases[e].local_nodes_for_primitive(lb.global_primitive_id(l), *state.mesh);
1015
1016 for (long n = 0; n < nodes.size(); ++n)
1017 {
1018 const assembler::AssemblyValues &v = vals.basis_values[nodes(n)];
1019 assert(v.global.size() == 1);
1020 for (int d = 0; d < actual_dim; d++)
1021 {
1022 double val = 0;
1023
1024 // j = j(x, grad u)
1025 if (j.depend_on_gradu())
1026 {
1027 for (int q = 0; q < weights.size(); ++q)
1028 val += dot(dj_dgradu.block(q, d * dim, 1, dim), v.grad_t_m.row(q));
1029 }
1030 // j = j(x, grad u)
1031 if (j.depend_on_gradu_local())
1032 {
1033 for (int q = 0; q < weights.size(); ++q)
1034 val += dot(dj_dgradu_local.block(q, d * dim, 1, dim), v.grad.row(q));
1035 }
1036 // j = j(x, u)
1037 if (j.depend_on_u())
1038 {
1039 for (int q = 0; q < weights.size(); ++q)
1040 val += dj_du(q, d) * v.val(q);
1041 }
1042 local_storage.vec(v.global[0].index * actual_dim + d) += val;
1043 }
1044 }
1045 }
1046 }
1047 }
1048 });
1049 for (const LocalThreadVecStorage &local_storage : storage)
1050 term += local_storage.vec;
1051 }
1052 else if (spatial_integral_type == SpatialIntegralType::VertexSum)
1053 {
1054 std::vector<bool> traversed(state.n_bases, false);
1055 IntegrableFunctional::ParameterType params;
1056 params.t = dt * cur_step + t0;
1057 params.step = cur_step;
1058 for (int e = 0; e < bases.size(); e++)
1059 {
1060 const auto &bs = bases[e];
1061 for (int i = 0; i < bs.bases.size(); i++)
1062 {
1063 const auto &b = bs.bases[i];
1064 assert(b.global().size() == 1);
1065 const auto &g = b.global()[0];
1066 if (traversed[g.index])
1067 continue;
1068
1069 const Eigen::MatrixXd lame_params = extract_lame_params(state.assembler->parameters(), e, params.t, Eigen::MatrixXd::Zero(1, dim) /*Not used*/, g.node);
1070
1071 params.node = g.index;
1072 params.elem = e;
1073 params.body_id = state.mesh->get_body_id(e);
1074 Eigen::MatrixXd val;
1075 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);
1076 term.block(g.index * actual_dim, 0, actual_dim, 1) += val.transpose();
1077 traversed[g.index] = true;
1078 }
1079 }
1080 }
1081 }
1082
1083 Eigen::VectorXd AdjointTools::map_primitive_to_node_order(const State &state, const Eigen::VectorXd &primitives)
1084 {
1085 int dim = state.mesh->dimension();
1086 assert(primitives.size() == (state.n_geom_bases * dim));
1087 Eigen::VectorXd nodes(primitives.size());
1088 auto map = state.primitive_to_node();
1089 for (int v = 0; v < state.n_geom_bases; ++v)
1090 nodes.segment(map[v] * dim, dim) = primitives.segment(v * dim, dim);
1091 return nodes;
1092 }
1093
1094 Eigen::VectorXd AdjointTools::map_node_to_primitive_order(const State &state, const Eigen::VectorXd &nodes)
1095 {
1096 int dim = state.mesh->dimension();
1097 assert(nodes.size() == (state.n_geom_bases * dim));
1098 Eigen::VectorXd primitives(nodes.size());
1099 auto map = state.node_to_primitive();
1100 for (int v = 0; v < state.n_geom_bases; ++v)
1101 primitives.segment(map[v] * dim, dim) = nodes.segment(v * dim, dim);
1102 return primitives;
1103 }
1104
1105 Eigen::MatrixXd AdjointTools::edge_normal_gradient(const Eigen::MatrixXd &V)
1106 {
1107 DiffScalarBase::setVariableCount(4);
1108 Eigen::Matrix<Diff, 4, 1> full_diff(4, 1);
1109 for (int i = 0; i < 2; i++)
1110 for (int j = 0; j < 2; j++)
1111 full_diff(i * 2 + j) = Diff(i * 2 + j, V(i, j));
1112 auto reduced_diff = edge_normal(full_diff);
1113
1114 Eigen::MatrixXd grad(2, 4);
1115 for (int i = 0; i < 2; ++i)
1116 grad.row(i) = reduced_diff[i].getGradient();
1117
1118 return grad;
1119 }
1120
1121 Eigen::MatrixXd AdjointTools::face_normal_gradient(const Eigen::MatrixXd &V)
1122 {
1123 DiffScalarBase::setVariableCount(9);
1124 Eigen::Matrix<Diff, 9, 1> full_diff(9, 1);
1125 for (int i = 0; i < 3; i++)
1126 for (int j = 0; j < 3; j++)
1127 full_diff(i * 3 + j) = Diff(i * 3 + j, V(i, j));
1128 auto reduced_diff = face_normal(full_diff);
1129
1130 Eigen::MatrixXd grad(3, 9);
1131 for (int i = 0; i < 3; ++i)
1132 grad.row(i) = reduced_diff[i].getGradient();
1133
1134 return grad;
1135 }
1136
1137 Eigen::MatrixXd AdjointTools::edge_velocity_divergence(const Eigen::MatrixXd &V)
1138 {
1139 return line_length_grad(V) / line_length<double>(V);
1140 }
1141
1142 Eigen::MatrixXd AdjointTools::face_velocity_divergence(const Eigen::MatrixXd &V)
1143 {
1144 return triangle_area_grad(V) / triangle_area<double>(V);
1145 }
1146} // namespace polyfem::solver
V
int V
Definition AdjointNLProblem.cpp:27
vec
Eigen::MatrixXd vec
Definition AdjointTools.cpp:71
val
double val
Definition AdjointTools.cpp:58
vals
assembler::ElementAssemblyValues vals
Definition AdjointTools.cpp:59
vec
Eigen::MatrixXd vec
Definition Assembler.cpp:72
val
double val
Definition Assembler.cpp:86
da
QuadratureVector da
Definition Assembler.cpp:23
vals
ElementAssemblyValues vals
Definition Assembler.cpp:22
gvals
assembler::ElementAssemblyValues gvals
Definition BodyForm.cpp:25
BodyForm.hpp
quadrature
Quadrature quadrature
Definition MassMatrixAssembler.cpp:30
polyfem::IntegrableFunctional::dj_du
void dj_du(const Eigen::MatrixXd &elastic_params, const Eigen::MatrixXd &local_pts, const Eigen::MatrixXd &pts, const Eigen::MatrixXd &u, const Eigen::MatrixXd &grad_u, const Eigen::MatrixXd &reference_normals, const assembler::ElementAssemblyValues &vals, ParameterType &params, Eigen::MatrixXd &val) const
Definition IntegrableFunctional.cpp:50
polyfem::IntegrableFunctional::dj_dgradu
void dj_dgradu(const Eigen::MatrixXd &elastic_params, const Eigen::MatrixXd &local_pts, const Eigen::MatrixXd &pts, const Eigen::MatrixXd &u, const Eigen::MatrixXd &grad_u, const Eigen::MatrixXd &reference_normals, const assembler::ElementAssemblyValues &vals, ParameterType &params, Eigen::MatrixXd &val) const
Definition IntegrableFunctional.cpp:56
polyfem::IntegrableFunctional::evaluate
void evaluate(const Eigen::MatrixXd &elastic_params, const Eigen::MatrixXd &local_pts, const Eigen::MatrixXd &pts, const Eigen::MatrixXd &u, const Eigen::MatrixXd &grad_u, const Eigen::MatrixXd &reference_normals, const assembler::ElementAssemblyValues &vals, ParameterType &params, Eigen::MatrixXd &val) const
Definition IntegrableFunctional.cpp:38
polyfem::IntegrableFunctional::dj_dx
void dj_dx(const Eigen::MatrixXd &elastic_params, const Eigen::MatrixXd &local_pts, const Eigen::MatrixXd &pts, const Eigen::MatrixXd &u, const Eigen::MatrixXd &grad_u, const Eigen::MatrixXd &reference_normals, const assembler::ElementAssemblyValues &vals, ParameterType &params, Eigen::MatrixXd &val) const
Definition IntegrableFunctional.cpp:44
polyfem::State
main class that contains the polyfem solver and all its state
Definition State.hpp:79
polyfem::State::n_bases
int n_bases
number of bases
Definition State.hpp:178
polyfem::State::ass_vals_cache
assembler::AssemblyValsCache ass_vals_cache
used to store assembly values for small problems
Definition State.hpp:196
polyfem::State::geom_bases
const std::vector< basis::ElementBases > & geom_bases() const
Get a constant reference to the geometry mapping bases.
Definition State.hpp:223
polyfem::State::assembler
std::shared_ptr< assembler::Assembler > assembler
assemblers
Definition State.hpp:155
polyfem::State::collision_mesh
ipc::CollisionMesh collision_mesh
IPC collision mesh.
Definition State.hpp:538
polyfem::State::mass_matrix_assembler
std::shared_ptr< assembler::Mass > mass_matrix_assembler
Definition State.hpp:157
polyfem::State::primitive_to_node
std::vector< int > primitive_to_node() const
Definition State.cpp:227
polyfem::State::mesh
std::unique_ptr< mesh::Mesh > mesh
current mesh, it can be a Mesh2D or Mesh3D
Definition State.hpp:466
polyfem::State::gbasis_nodes_to_basis_nodes
StiffnessMatrix gbasis_nodes_to_basis_nodes
Definition State.hpp:719
polyfem::State::problem
std::shared_ptr< assembler::Problem > problem
current problem, it contains rhs and bc
Definition State.hpp:168
polyfem::State::node_to_primitive
std::vector< int > node_to_primitive() const
Definition State.cpp:234
polyfem::State::args
json args
main input arguments containing all defaults
Definition State.hpp:101
polyfem::State::diff_cached
solver::DiffCache diff_cached
Definition State.hpp:669
polyfem::State::bases
std::vector< basis::ElementBases > bases
FE bases, the size is #elements.
Definition State.hpp:171
polyfem::State::n_geom_bases
int n_geom_bases
number of geometric bases
Definition State.hpp:182
polyfem::State::total_local_boundary
std::vector< mesh::LocalBoundary > total_local_boundary
mapping from elements to nodes for all mesh
Definition State.hpp:431
polyfem::State::ndof
int ndof() const
Definition State.hpp:673
polyfem::State::mass_ass_vals_cache
assembler::AssemblyValsCache mass_ass_vals_cache
Definition State.hpp:197
polyfem::State::boundary_nodes
std::vector< int > boundary_nodes
list of boundary nodes
Definition State.hpp:427
polyfem::State::solve_data
solver::SolveData solve_data
timedependent stuff cached
Definition State.hpp:324
polyfem::State::is_contact_enabled
bool is_contact_enabled() const
does the simulation has contact
Definition State.hpp:574
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:33
polyfem::assembler::AssemblyValues
stores per local bases evaluations
Definition AssemblyValues.hpp:14
polyfem::assembler::AssemblyValues::global
std::vector< basis::Local2Global > global
Definition AssemblyValues.hpp:20
polyfem::assembler::ElementAssemblyValues
stores per element basis values at given quadrature points and geometric mapping
Definition ElementAssemblyValues.hpp:14
polyfem::assembler::ElementAssemblyValues::compute
void compute(const int el_index, const bool is_volume, const Eigen::MatrixXd &pts, const basis::ElementBases &basis, const basis::ElementBases &gbasis)
computes the per element values at the local (ref el) points (pts) sets basis_values,...
Definition ElementAssemblyValues.cpp:151
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:983
polyfem::solver::DiffCache::collision_set
const ipc::Collisions & collision_set(int step) const
Definition DiffCache.hpp:143
polyfem::solver::DiffCache::v
Eigen::VectorXd v(int step) const
Definition DiffCache.hpp:119
polyfem::solver::DiffCache::u
Eigen::VectorXd u(int step) const
Definition DiffCache.hpp:112
polyfem::solver::DiffCache::friction_collision_set
const ipc::FrictionCollisions & friction_collision_set(int step) const
Definition DiffCache.hpp:150
polyfem::solver::SolveData::friction_form
std::shared_ptr< solver::FrictionForm > friction_form
Definition SolveData.hpp:150
polyfem::solver::SolveData::inertia_form
std::shared_ptr< solver::InertiaForm > inertia_form
Definition SolveData.hpp:151
polyfem::solver::SolveData::pressure_form
std::shared_ptr< solver::PressureForm > pressure_form
Definition SolveData.hpp:152
polyfem::solver::SolveData::body_form
std::shared_ptr< solver::BodyForm > body_form
Definition SolveData.hpp:146
polyfem::solver::SolveData::contact_form
std::shared_ptr< solver::ContactForm > contact_form
Definition SolveData.hpp:147
polyfem::solver::SolveData::damping_form
std::shared_ptr< solver::ElasticForm > damping_form
Definition SolveData.hpp:148
polyfem::solver::SolveData::elastic_form
std::shared_ptr< solver::ElasticForm > elastic_form
Definition SolveData.hpp:149
generate_mooney_rivlin.grad
grad
Definition generate_mooney_rivlin.py:51
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:357
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:177
polyfem::utils::vector2matrix
void vector2matrix(const Eigen::VectorXd &vec, Eigen::MatrixXd &mat)
Definition MatrixUtils.cpp:69
polyfem::utils::get_local_thread_storage
auto & get_local_thread_storage(Storages &storage, int thread_id)
polyfem::utils::create_thread_storage
auto create_thread_storage(const LocalStorage &initial_local_storage)
polyfem::utils::triangle_area
double triangle_area(const Eigen::MatrixXd V)
Compute the signed area of a triangle defined by three points.
Definition GeometryUtils.cpp:30
polyfem::utils::maybe_parallel_for
void maybe_parallel_for(int size, const std::function< void(int, int, int)> &partial_for)
polyfem::QuadratureVector
Eigen::Matrix< double, Eigen::Dynamic, 1, 0, MAX_QUAD_POINTS, 1 > QuadratureVector
Definition Types.hpp:15
polyfem::log_and_throw_adjoint_error
void log_and_throw_adjoint_error(const std::string &msg)
Definition Logger.cpp:77
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