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