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