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