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