Loading [MathJax]/extensions/tex2jax.js
PolyFEM
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
AMIPSEnergy.cpp
Go to the documentation of this file.
1#include "AMIPSEnergy.hpp"
2
3#include <polyfem/utils/Logger.hpp>
4
5namespace polyfem::assembler
6{
7 namespace
8 {
9
10 bool delta(int i, int j)
11 {
12 return (i == j) ? true : false;
13 }
14
15 template <class T, int p = 3>
16 class barrier
17 {
18 constexpr static double C = 1e6;
19
20 public:
21 static_assert(p % 2 == 1);
22 static T value(T J)
23 {
24 if (J >= 0.5)
25 return T(0.);
26 const T tmp1 = 2 * J - 1;
27 const T tmp2 = pow(tmp1, p);
28 return C * (1 / (tmp2 + 1) - 1);
29 }
30
31 static T first_derivatives(T J)
32 {
33 if (J >= 0.5)
34 return T(0);
35 const T tmp1 = 2 * J - 1;
36 const T tmp2 = pow(tmp1, p);
37 return C * -2 * p * tmp2 / tmp1 / pow(1 + tmp2, 2);
38 }
39
40 static T second_derivatives(T J)
41 {
42 if (J >= 0.5)
43 return T(0);
44 const T tmp1 = 2 * J - 1;
45 const T tmp2 = pow(tmp1, p);
46 return C * 4 * p * tmp2 / pow(tmp1, 2) * ((1 - p) + (1 + p) * tmp2) / pow(1 + tmp2, 3);
47 }
48 };
49
50 template <int dim>
51 Eigen::Matrix<double, dim, dim> hat(const Eigen::Matrix<double, dim, 1> &x)
52 {
53
54 Eigen::Matrix<double, dim, dim> prod;
55 prod.setZero();
56
57 prod(0, 1) = -x(2);
58 prod(0, 2) = x(1);
59 prod(1, 0) = x(2);
60 prod(1, 2) = -x(0);
61 prod(2, 0) = -x(1);
62 prod(2, 1) = x(0);
63
64 return prod;
65 }
66
67 template <int dim>
68 Eigen::Matrix<double, dim, 1> cross(const Eigen::Matrix<double, dim, 1> &x, const Eigen::Matrix<double, dim, 1> &y)
69 {
70
71 Eigen::Matrix<double, dim, 1> z;
72 z.setZero();
73
74 z(0) = x(1) * y(2) - x(2) * y(1);
75 z(1) = x(2) * y(0) - x(0) * y(2);
76 z(2) = x(0) * y(1) - x(1) * y(0);
77
78 return z;
79 }
80
81 template <int dimt, class T>
82 Eigen::Matrix<T, dimt, dimt> get_standard(const int dim, const bool use_rest_pose)
83 {
84 Eigen::Matrix<double, dimt, dimt> standard(dim, dim);
85 if (use_rest_pose)
86 {
87 standard.setIdentity();
88 }
89 else
90 {
91 if (dim == 2)
92 standard << 1, 0,
93 0.5, std::sqrt(3) / 2;
94 else
95 standard << 1, 0, 0,
96 0.5, std::sqrt(3) / 2., 0,
97 0.5, 0.5 / std::sqrt(3), std::sqrt(3) / 2.;
98 standard = standard.inverse().transpose().eval();
99 }
100
101 Eigen::Matrix<T, dimt, dimt> res(dim, dim);
102 for (int i = 0; i < dim; ++i)
103 {
104 for (int j = 0; j < dim; ++j)
105 {
106 res(i, j) = T(standard(i, j));
107 }
108 }
109
110 return res;
111 }
112 } // namespace
113 double AMIPSEnergy::compute_energy(const NonLinearAssemblerData &data) const
114 {
115 return compute_energy_aux<double>(data);
116 }
117
118 Eigen::VectorXd AMIPSEnergy::assemble_gradient(const NonLinearAssemblerData &data) const
119 {
120 const int n_bases = data.vals.basis_values.size();
121 return polyfem::gradient_from_energy(
122 size(), n_bases, data,
123 [&](const NonLinearAssemblerData &data) { return compute_energy_aux<DScalar1<double, Eigen::Matrix<double, 6, 1>>>(data); },
124 [&](const NonLinearAssemblerData &data) { return compute_energy_aux<DScalar1<double, Eigen::Matrix<double, 8, 1>>>(data); },
125 [&](const NonLinearAssemblerData &data) { return compute_energy_aux<DScalar1<double, Eigen::Matrix<double, 12, 1>>>(data); },
126 [&](const NonLinearAssemblerData &data) { return compute_energy_aux<DScalar1<double, Eigen::Matrix<double, 18, 1>>>(data); },
127 [&](const NonLinearAssemblerData &data) { return compute_energy_aux<DScalar1<double, Eigen::Matrix<double, 24, 1>>>(data); },
128 [&](const NonLinearAssemblerData &data) { return compute_energy_aux<DScalar1<double, Eigen::Matrix<double, 30, 1>>>(data); },
129 [&](const NonLinearAssemblerData &data) { return compute_energy_aux<DScalar1<double, Eigen::Matrix<double, 60, 1>>>(data); },
130 [&](const NonLinearAssemblerData &data) { return compute_energy_aux<DScalar1<double, Eigen::Matrix<double, 81, 1>>>(data); },
131 [&](const NonLinearAssemblerData &data) { return compute_energy_aux<DScalar1<double, Eigen::Matrix<double, Eigen::Dynamic, 1, 0, SMALL_N, 1>>>(data); },
132 [&](const NonLinearAssemblerData &data) { return compute_energy_aux<DScalar1<double, Eigen::Matrix<double, Eigen::Dynamic, 1, 0, BIG_N, 1>>>(data); },
133 [&](const NonLinearAssemblerData &data) { return compute_energy_aux<DScalar1<double, Eigen::VectorXd>>(data); });
134
135 // if (size() == 2)
136 // {
137 // switch (data.vals.basis_values.size())
138 // {
139 // case 3:
140 // {
141 // gradient.resize(6);
142 // compute_energy_aux_gradient_fast<3, 2>(data, gradient);
143 // break;
144 // }
145 // case 6:
146 // {
147 // gradient.resize(12);
148 // compute_energy_aux_gradient_fast<6, 2>(data, gradient);
149 // break;
150 // }
151 // case 10:
152 // {
153 // gradient.resize(20);
154 // compute_energy_aux_gradient_fast<10, 2>(data, gradient);
155 // break;
156 // }
157 // default:
158 // {
159 // gradient.resize(data.vals.basis_values.size() * 2);
160 // compute_energy_aux_gradient_fast<Eigen::Dynamic, 2>(data, gradient);
161 // break;
162 // }
163 // }
164 // }
165 // else // if (size() == 3)
166 // {
167 // assert(size() == 3);
168 // switch (data.vals.basis_values.size())
169 // {
170 // case 4:
171 // {
172 // gradient.resize(12);
173 // compute_energy_aux_gradient_fast<4, 3>(data, gradient);
174 // break;
175 // }
176 // case 10:
177 // {
178 // gradient.resize(30);
179 // compute_energy_aux_gradient_fast<10, 3>(data, gradient);
180 // break;
181 // }
182 // case 20:
183 // {
184 // gradient.resize(60);
185 // compute_energy_aux_gradient_fast<20, 3>(data, gradient);
186 // break;
187 // }
188 // default:
189 // {
190 // gradient.resize(data.vals.basis_values.size() * 3);
191 // compute_energy_aux_gradient_fast<Eigen::Dynamic, 3>(data, gradient);
192 // break;
193 // }
194 // }
195 // }
196 }
197
198 // Compute ∫ tr(FᵀF) / J^(1+2/dim) dxdydz
199 template <typename T>
200 T AMIPSEnergy::compute_energy_aux(const NonLinearAssemblerData &data) const
201 {
202 typedef Eigen::Matrix<T, Eigen::Dynamic, 1> AutoDiffVect;
203 typedef Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic, 0, 3, 3> AutoDiffGradMat;
204 typedef Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, 0, 3, 3> DoubleGradMat;
205
206 double power = -1;
207 if (use_rest_pose_)
208 power = size() == 2 ? 1. : (2. / 3.);
209 else
210 power = size() == 2 ? 2. : 5. / 3.;
211
212 AutoDiffVect local_disp;
213 get_local_disp(data, size(), local_disp);
214
215 AutoDiffGradMat standard;
216
217 if (size() == 2)
218 standard = get_standard<2, T>(size(), use_rest_pose_);
219 else
220 standard = get_standard<3, T>(size(), use_rest_pose_);
221
222 AutoDiffGradMat def_grad(size(), size());
223
224 T energy = T(0.0);
225
226 const int n_pts = data.da.size();
227 for (long p = 0; p < n_pts; ++p)
228 {
229 compute_disp_grad_at_quad(data, local_disp, p, size(), def_grad);
230
231 // Id + grad d
232 for (int d = 0; d < size(); ++d)
233 def_grad(d, d) += T(1);
234
235 if (!use_rest_pose_)
236 {
237 DoubleGradMat tmp_jac_it = data.vals.jac_it[p];
238 tmp_jac_it = tmp_jac_it.inverse();
239
240 AutoDiffGradMat jac_it(size(), size());
241 for (long k = 0; k < jac_it.size(); ++k)
242 jac_it(k) = T(tmp_jac_it(k));
243
244 def_grad *= jac_it;
245 def_grad = def_grad * standard;
246 }
247
248 const T det = polyfem::utils::determinant(def_grad);
249 if (det <= 0)
250 {
251 energy = std::nan("");
252 break;
253 }
254
255 const T powJ = pow(det, power);
256 const T val = (def_grad.transpose() * def_grad).trace() / powJ; //+ barrier<T>::value(det);
257
258 energy += val * data.da(p);
259 }
260 return energy;
261 }
262
263 template <int n_basis, int dim>
264 void AMIPSEnergy::compute_energy_aux_gradient_fast(const NonLinearAssemblerData &data, Eigen::Matrix<double, Eigen::Dynamic, 1> &G_flattened) const
265 {
266 assert(data.x.cols() == 1);
267
268 const int n_pts = data.da.size();
269
270 Eigen::Matrix<double, n_basis, dim> local_disp(data.vals.basis_values.size(), size());
271 local_disp.setZero();
272 for (size_t i = 0; i < data.vals.basis_values.size(); ++i)
273 {
274 const auto &bs = data.vals.basis_values[i];
275 for (size_t ii = 0; ii < bs.global.size(); ++ii)
276 {
277 for (int d = 0; d < size(); ++d)
278 {
279 local_disp(i, d) += bs.global[ii].val * data.x(bs.global[ii].index * size() + d);
280 }
281 }
282 }
283
284 Eigen::Matrix<double, dim, dim> def_grad(size(), size());
285
286 Eigen::Matrix<double, n_basis, dim> G(data.vals.basis_values.size(), size());
287 G.setZero();
288
289 const Eigen::Matrix<double, dim, dim> standard = get_standard<dim, double>(size(), use_rest_pose_);
290
291 for (long p = 0; p < n_pts; ++p)
292 {
293 Eigen::Matrix<double, n_basis, dim> grad(data.vals.basis_values.size(), size());
294
295 for (size_t i = 0; i < data.vals.basis_values.size(); ++i)
296 {
297 grad.row(i) = data.vals.basis_values[i].grad.row(p);
298 }
299
300 Eigen::Matrix<double, dim, dim> jac_it = data.vals.jac_it[p];
301
302 // Id + grad d
303 def_grad = (local_disp.transpose() * grad + jac_it.inverse()) * standard;
304
305 double J = def_grad.determinant();
306 if (J <= 0)
307 J = std::nan("");
308
309 Eigen::Matrix<double, dim, dim> delJ_delF(size(), size());
310 delJ_delF.setZero();
311
312 if (dim == 2)
313 {
314
315 delJ_delF(0, 0) = def_grad(1, 1);
316 delJ_delF(0, 1) = -def_grad(1, 0);
317 delJ_delF(1, 0) = -def_grad(0, 1);
318 delJ_delF(1, 1) = def_grad(0, 0);
319 }
320
321 else if (dim == 3)
322 {
323 Eigen::Matrix<double, dim, 1> u(def_grad.rows());
324 Eigen::Matrix<double, dim, 1> v(def_grad.rows());
325 Eigen::Matrix<double, dim, 1> w(def_grad.rows());
326
327 u = def_grad.col(0);
328 v = def_grad.col(1);
329 w = def_grad.col(2);
330
331 delJ_delF.col(0) = cross<dim>(v, w);
332 delJ_delF.col(1) = cross<dim>(w, u);
333 delJ_delF.col(2) = cross<dim>(u, v);
334 }
335
336 const double m = (dim == 2) ? 2. : 5. / 3.;
337 const double powJ = pow(J, m);
338 Eigen::Matrix<double, dim, dim> gradient_temp = (2 * def_grad - (def_grad.squaredNorm() * m) / J * delJ_delF) / powJ;
339 Eigen::Matrix<double, n_basis, dim> gradient = grad * standard * gradient_temp.transpose();
340
341 G.noalias() += gradient * data.da(p);
342 }
343
344 Eigen::Matrix<double, dim, n_basis> G_T = G.transpose();
345
346 constexpr int N = (n_basis == Eigen::Dynamic) ? Eigen::Dynamic : n_basis * dim;
347 Eigen::Matrix<double, N, 1> temp(Eigen::Map<Eigen::Matrix<double, N, 1>>(G_T.data(), G_T.size()));
348 G_flattened = temp;
349 }
350
351 Eigen::MatrixXd AMIPSEnergy::assemble_hessian(const NonLinearAssemblerData &data) const
352 {
353 const int n_bases = data.vals.basis_values.size();
354 return polyfem::hessian_from_energy(
355 size(), n_bases, data,
356 [&](const NonLinearAssemblerData &data) { return compute_energy_aux<DScalar2<double, Eigen::Matrix<double, 6, 1>, Eigen::Matrix<double, 6, 6>>>(data); },
357 [&](const NonLinearAssemblerData &data) { return compute_energy_aux<DScalar2<double, Eigen::Matrix<double, 8, 1>, Eigen::Matrix<double, 8, 8>>>(data); },
358 [&](const NonLinearAssemblerData &data) { return compute_energy_aux<DScalar2<double, Eigen::Matrix<double, 12, 1>, Eigen::Matrix<double, 12, 12>>>(data); },
359 [&](const NonLinearAssemblerData &data) { return compute_energy_aux<DScalar2<double, Eigen::Matrix<double, 18, 1>, Eigen::Matrix<double, 18, 18>>>(data); },
360 [&](const NonLinearAssemblerData &data) { return compute_energy_aux<DScalar2<double, Eigen::Matrix<double, 24, 1>, Eigen::Matrix<double, 24, 24>>>(data); },
361 [&](const NonLinearAssemblerData &data) { return compute_energy_aux<DScalar2<double, Eigen::Matrix<double, 30, 1>, Eigen::Matrix<double, 30, 30>>>(data); },
362 [&](const NonLinearAssemblerData &data) { return compute_energy_aux<DScalar2<double, Eigen::Matrix<double, 60, 1>, Eigen::Matrix<double, 60, 60>>>(data); },
363 [&](const NonLinearAssemblerData &data) { return compute_energy_aux<DScalar2<double, Eigen::Matrix<double, 81, 1>, Eigen::Matrix<double, 81, 81>>>(data); },
364 [&](const NonLinearAssemblerData &data) { return compute_energy_aux<DScalar2<double, Eigen::Matrix<double, Eigen::Dynamic, 1, 0, SMALL_N, 1>, Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, 0, SMALL_N, SMALL_N>>>(data); },
365 [&](const NonLinearAssemblerData &data) { return compute_energy_aux<DScalar2<double, Eigen::VectorXd, Eigen::MatrixXd>>(data); });
366
367 // Eigen::MatrixXd hessian;
368
369 // if (size() == 2)
370 // {
371 // switch (data.vals.basis_values.size())
372 // {
373 // case 3:
374 // {
375 // hessian.resize(6, 6);
376 // hessian.setZero();
377 // compute_energy_hessian_aux_fast<3, 2>(data, hessian);
378 // break;
379 // }
380 // case 6:
381 // {
382 // hessian.resize(12, 12);
383 // hessian.setZero();
384 // compute_energy_hessian_aux_fast<6, 2>(data, hessian);
385 // break;
386 // }
387 // case 10:
388 // {
389 // hessian.resize(20, 20);
390 // hessian.setZero();
391 // compute_energy_hessian_aux_fast<10, 2>(data, hessian);
392 // break;
393 // }
394 // default:
395 // {
396 // hessian.resize(data.vals.basis_values.size() * 2, data.vals.basis_values.size() * 2);
397 // hessian.setZero();
398 // compute_energy_hessian_aux_fast<Eigen::Dynamic, 2>(data, hessian);
399 // break;
400 // }
401 // }
402 // }
403 // else // if (size() == 3)
404 // {
405 // assert(size() == 3);
406 // switch (data.vals.basis_values.size())
407 // {
408 // case 4:
409 // {
410 // hessian.resize(12, 12);
411 // hessian.setZero();
412 // compute_energy_hessian_aux_fast<4, 3>(data, hessian);
413 // break;
414 // }
415 // case 10:
416 // {
417 // hessian.resize(30, 30);
418 // hessian.setZero();
419 // compute_energy_hessian_aux_fast<10, 3>(data, hessian);
420 // break;
421 // }
422 // case 20:
423 // {
424 // hessian.resize(60, 60);
425 // hessian.setZero();
426 // compute_energy_hessian_aux_fast<20, 3>(data, hessian);
427 // break;
428 // }
429 // default:
430 // {
431 // hessian.resize(data.vals.basis_values.size() * 3, data.vals.basis_values.size() * 3);
432 // hessian.setZero();
433 // compute_energy_hessian_aux_fast<Eigen::Dynamic, 3>(data, hessian);
434 // break;
435 // }
436 // }
437 // }
438
439 // return hessian;
440 }
441
442 template <int n_basis, int dim>
443 void AMIPSEnergy::compute_energy_hessian_aux_fast(const NonLinearAssemblerData &data, Eigen::MatrixXd &H) const
444 {
445 assert(data.x.cols() == 1);
446
447 constexpr int N = (n_basis == Eigen::Dynamic) ? Eigen::Dynamic : n_basis * dim;
448 const int n_pts = data.da.size();
449
450 Eigen::Matrix<double, n_basis, dim> local_disp(data.vals.basis_values.size(), size());
451 local_disp.setZero();
452 for (size_t i = 0; i < data.vals.basis_values.size(); ++i)
453 {
454 const auto &bs = data.vals.basis_values[i];
455 for (size_t ii = 0; ii < bs.global.size(); ++ii)
456 {
457 for (int d = 0; d < size(); ++d)
458 {
459 local_disp(i, d) += bs.global[ii].val * data.x(bs.global[ii].index * size() + d);
460 }
461 }
462 }
463
464 Eigen::Matrix<double, dim, dim> def_grad(size(), size());
465
466 const Eigen::Matrix<double, dim, dim> standard = get_standard<dim, double>(size(), use_rest_pose_);
467
468 for (long p = 0; p < n_pts; ++p)
469 {
470 Eigen::Matrix<double, n_basis, dim> grad(data.vals.basis_values.size(), size());
471
472 for (size_t i = 0; i < data.vals.basis_values.size(); ++i)
473 {
474 grad.row(i) = data.vals.basis_values[i].grad.row(p);
475 }
476
477 Eigen::Matrix<double, dim, dim> jac_it = data.vals.jac_it[p];
478
479 // Id + grad d
480 if (use_rest_pose_)
481 {
482 def_grad = local_disp.transpose() * grad + Eigen::Matrix<double, dim, dim>::Identity(size(), size());
483 }
484 else
485 {
486 def_grad = (local_disp.transpose() * grad + jac_it.inverse()) * standard;
487 }
488
489 Eigen::Matrix<double, dim * dim, dim * dim> hessian_temp;
490 {
491 typedef DScalar2<double, Eigen::Matrix<double, dim * dim, 1>, Eigen::Matrix<double, dim * dim, dim * dim>> Diff;
492 DiffScalarBase::setVariableCount(dim * dim);
493
494 Eigen::Matrix<Diff, dim, dim> def_grad_ad(dim, dim);
495 for (int i = 0; i < def_grad.size(); i++)
496 def_grad_ad(i) = Diff(i, def_grad(i));
497
498 Diff val = pow(utils::determinant(def_grad_ad), -1. - 2. / dim) * def_grad_ad.squaredNorm();
499
500 hessian_temp = val.getHessian();
501 }
502
503 // const double J = def_grad.determinant();
504 // const double TrFFt = def_grad.squaredNorm();
505
506 // Eigen::Matrix<double, dim, dim> delJ_delF(size(), size());
507 // delJ_delF.setZero();
508 // Eigen::Matrix<double, dim * dim, dim * dim> del2J_delF2(size() * size(), size() * size());
509 // del2J_delF2.setZero();
510
511 // if (dim == 2)
512 // {
513 // delJ_delF(0, 0) = def_grad(1, 1);
514 // delJ_delF(0, 1) = -def_grad(1, 0);
515 // delJ_delF(1, 0) = -def_grad(0, 1);
516 // delJ_delF(1, 1) = def_grad(0, 0);
517
518 // del2J_delF2(0, 3) = 1;
519 // del2J_delF2(1, 2) = -1;
520 // del2J_delF2(2, 1) = -1;
521 // del2J_delF2(3, 0) = 1;
522 // }
523 // else if (size() == 3)
524 // {
525 // Eigen::Matrix<double, dim, 1> u(def_grad.rows());
526 // Eigen::Matrix<double, dim, 1> v(def_grad.rows());
527 // Eigen::Matrix<double, dim, 1> w(def_grad.rows());
528
529 // u = def_grad.col(0);
530 // v = def_grad.col(1);
531 // w = def_grad.col(2);
532
533 // delJ_delF.col(0) = cross<dim>(v, w);
534 // delJ_delF.col(1) = cross<dim>(w, u);
535 // delJ_delF.col(2) = cross<dim>(u, v);
536
537 // del2J_delF2.template block<dim, dim>(0, 6) = hat<dim>(v);
538 // del2J_delF2.template block<dim, dim>(6, 0) = -hat<dim>(v);
539 // del2J_delF2.template block<dim, dim>(0, 3) = -hat<dim>(w);
540 // del2J_delF2.template block<dim, dim>(3, 0) = hat<dim>(w);
541 // del2J_delF2.template block<dim, dim>(3, 6) = -hat<dim>(u);
542 // del2J_delF2.template block<dim, dim>(6, 3) = hat<dim>(u);
543 // }
544
545 // Eigen::Matrix<double, dim * dim, dim * dim> id = Eigen::Matrix<double, dim * dim, dim * dim>::Identity(size() * size(), size() * size());
546
547 // Eigen::Matrix<double, dim * dim, 1> g_j = Eigen::Map<const Eigen::Matrix<double, dim * dim, 1>>(delJ_delF.data(), delJ_delF.size());
548 // Eigen::Matrix<double, dim * dim, 1> F_flattened = Eigen::Map<const Eigen::Matrix<double, dim * dim, 1>>(def_grad.data(), def_grad.size());
549
550 // const double tmp = TrFFt / J / dim;
551 // const double Jpow = pow(J, 2. / dim);
552 // Eigen::Matrix<double, dim * dim, dim * dim> hessian_temp = -4. / dim / (Jpow * J) * (F_flattened - tmp * g_j) * g_j.transpose() +
553 // 2. / Jpow * (id - ((2 / dim * F_flattened - tmp * g_j) / J * g_j.transpose() + tmp * del2J_delF2));
554
555 Eigen::Matrix<double, dim * dim, N> delF_delU_tensor(jac_it.size(), grad.size());
556
557 for (size_t i = 0; i < local_disp.rows(); ++i)
558 {
559 for (size_t j = 0; j < local_disp.cols(); ++j)
560 {
561 Eigen::Matrix<double, dim, dim> temp(size(), size());
562 temp.setZero();
563 temp.row(j) = grad.row(i);
564 temp = temp * standard;
565 Eigen::Matrix<double, dim * dim, 1> temp_flattened(Eigen::Map<Eigen::Matrix<double, dim * dim, 1>>(temp.data(), temp.size()));
566 delF_delU_tensor.col(i * size() + j) = temp_flattened;
567 }
568 }
569
570 Eigen::Matrix<double, N, N> hessian = delF_delU_tensor.transpose() * hessian_temp * delF_delU_tensor;
571
572 H += hessian * data.da(p);
573 }
574 }
575
576 void AMIPSEnergy::assign_stress_tensor(const OutputData &data,
577 const int all_size,
578 const ElasticityTensorType &type,
579 Eigen::MatrixXd &all,
580 const std::function<Eigen::MatrixXd(const Eigen::MatrixXd &)> &fun) const
581 {
582 const auto &displacement = data.fun;
583 const auto &local_pts = data.local_pts;
584 const auto &bs = data.bs;
585 const auto &gbs = data.gbs;
586 const auto el_id = data.el_id;
587 const auto t = data.t;
588
589 Eigen::MatrixXd displacement_grad(size(), size());
590
591 assert(displacement.cols() == 1);
592
593 all.setZero(local_pts.rows(), all_size);
594 }
595
596 void AMIPSEnergy::add_multimaterial(const int index, const json &params, const Units &units)
597 {
598 assert(size() == 2 || size() == 3);
599
600 if (params.contains("use_rest_pose"))
601 {
602 use_rest_pose_ = params["use_rest_pose"].get<bool>();
603 }
604 }
605
610
611 std::map<std::string, Assembler::ParamFunc> AMIPSEnergyAutodiff::parameters() const
612 {
613 return std::map<std::string, ParamFunc>();
614 }
615
616 void AMIPSEnergyAutodiff::add_multimaterial(const int index, const json &params, const Units &)
617 {
618 if (params.contains("canonical_transformation"))
619 {
620 canonical_transformation_.reserve(params["canonical_transformation"].size());
621 for (int i = 0; i < params["canonical_transformation"].size(); ++i)
622 {
623 Eigen::MatrixXd transform_matrix(size(), size());
624 for (int j = 0; j < size(); ++j)
625 for (int k = 0; k < size(); ++k)
626 transform_matrix(j, k) = params["canonical_transformation"][i][j][k];
627 canonical_transformation_.push_back(transform_matrix);
628 }
629 }
630 }
631
632} // namespace polyfem::assembler
val
double val
Definition Assembler.cpp:86
y
int y
Definition SplineBasis3d.cpp:55
z
int z
Definition SplineBasis3d.cpp:55
x
int x
Definition SplineBasis3d.cpp:55
polyfem::assembler::AMIPSEnergyAutodiff::parameters
std::map< std::string, ParamFunc > parameters() const override
Definition AMIPSEnergy.cpp:611
polyfem::assembler::AMIPSEnergyAutodiff::canonical_transformation_
std::vector< Eigen::MatrixXd > canonical_transformation_
Definition AMIPSEnergy.hpp:98
polyfem::assembler::AMIPSEnergyAutodiff::add_multimaterial
void add_multimaterial(const int index, const json &params, const Units &units) override
Definition AMIPSEnergy.cpp:616
polyfem::assembler::AMIPSEnergy::compute_energy
double compute_energy(const NonLinearAssemblerData &data) const override
Definition AMIPSEnergy.cpp:113
polyfem::assembler::AMIPSEnergy::compute_energy_aux
T compute_energy_aux(const NonLinearAssemblerData &data) const
Definition AMIPSEnergy.cpp:200
polyfem::assembler::AMIPSEnergy::assign_stress_tensor
void assign_stress_tensor(const OutputData &data, const int all_size, const ElasticityTensorType &type, Eigen::MatrixXd &all, const std::function< Eigen::MatrixXd(const Eigen::MatrixXd &)> &fun) const override
Definition AMIPSEnergy.cpp:576
polyfem::assembler::AMIPSEnergy::compute_energy_hessian_aux_fast
void compute_energy_hessian_aux_fast(const NonLinearAssemblerData &data, Eigen::MatrixXd &H) const
Definition AMIPSEnergy.cpp:443
polyfem::assembler::AMIPSEnergy::add_multimaterial
void add_multimaterial(const int index, const json &params, const Units &units) override
Definition AMIPSEnergy.cpp:596
polyfem::assembler::AMIPSEnergy::assemble_hessian
Eigen::MatrixXd assemble_hessian(const NonLinearAssemblerData &data) const override
Definition AMIPSEnergy.cpp:351
polyfem::assembler::AMIPSEnergy::compute_energy_aux_gradient_fast
void compute_energy_aux_gradient_fast(const NonLinearAssemblerData &data, Eigen::VectorXd &G_flattened) const
polyfem::assembler::AMIPSEnergy::assemble_gradient
Eigen::VectorXd assemble_gradient(const NonLinearAssemblerData &data) const override
Definition AMIPSEnergy.cpp:118
polyfem::assembler::ElementAssemblyValues::jac_it
std::vector< Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic, 0, 3, 3 > > jac_it
Definition ElementAssemblyValues.hpp:22
polyfem::assembler::OutputData::bs
const basis::ElementBases & bs
Definition AssemblerData.hpp:119
polyfem::assembler::OutputData::fun
const Eigen::MatrixXd & fun
Definition AssemblerData.hpp:122
polyfem::assembler::OutputData::local_pts
const Eigen::MatrixXd & local_pts
Definition AssemblerData.hpp:121
polyfem::assembler::OutputData::gbs
const basis::ElementBases & gbs
Definition AssemblerData.hpp:120
polyfem::assembler
Used for test only.
Definition AMIPSEnergy.cpp:6
polyfem::assembler::AutoDiffGradMat
Eigen::Matrix< T, Eigen::Dynamic, Eigen::Dynamic, 0, 3, 3 > AutoDiffGradMat
Definition MooneyRivlin3ParamSymbolic.hpp:12
polyfem::assembler::hat
Eigen::Matrix< double, dim, dim > hat(const Eigen::Matrix< double, dim, 1 > &x)
Definition NeoHookeanElasticity.cpp:422
polyfem::assembler::cross
Eigen::Matrix< double, dim, 1 > cross(const Eigen::Matrix< double, dim, 1 > &x, const Eigen::Matrix< double, dim, 1 > &y)
Definition NeoHookeanElasticity.cpp:439
polyfem::assembler::AutoDiffVect
Eigen::Matrix< T, Eigen::Dynamic, 1 > AutoDiffVect
Definition MooneyRivlin3ParamSymbolic.hpp:14
polyfem::solver::Diff
DScalar1< double, Eigen::Matrix< double, Eigen::Dynamic, 1 > > Diff
Definition SurfaceTractionForms.cpp:235
polyfem::utils::determinant
T determinant(const Eigen::Matrix< T, rows, cols, option, maxRow, maxCol > &mat)
Definition MatrixUtils.hpp:22
polyfem::compute_disp_grad_at_quad
void compute_disp_grad_at_quad(const assembler::NonLinearAssemblerData &data, const AutoDiffVect &local_disp, const int p, const int size, AutoDiffGradMat &def_grad)
Definition ElasticityUtils.hpp:106
polyfem::hessian_from_energy
Eigen::MatrixXd hessian_from_energy(const int size, const int n_bases, const assembler::NonLinearAssemblerData &data, const std::function< DScalar2< double, Eigen::Matrix< double, 6, 1 >, Eigen::Matrix< double, 6, 6 > >(const assembler::NonLinearAssemblerData &)> &fun6, const std::function< DScalar2< double, Eigen::Matrix< double, 8, 1 >, Eigen::Matrix< double, 8, 8 > >(const assembler::NonLinearAssemblerData &)> &fun8, const std::function< DScalar2< double, Eigen::Matrix< double, 12, 1 >, Eigen::Matrix< double, 12, 12 > >(const assembler::NonLinearAssemblerData &)> &fun12, const std::function< DScalar2< double, Eigen::Matrix< double, 18, 1 >, Eigen::Matrix< double, 18, 18 > >(const assembler::NonLinearAssemblerData &)> &fun18, const std::function< DScalar2< double, Eigen::Matrix< double, 24, 1 >, Eigen::Matrix< double, 24, 24 > >(const assembler::NonLinearAssemblerData &)> &fun24, const std::function< DScalar2< double, Eigen::Matrix< double, 30, 1 >, Eigen::Matrix< double, 30, 30 > >(const assembler::NonLinearAssemblerData &)> &fun30, const std::function< DScalar2< double, Eigen::Matrix< double, 60, 1 >, Eigen::Matrix< double, 60, 60 > >(const assembler::NonLinearAssemblerData &)> &fun60, const std::function< DScalar2< double, Eigen::Matrix< double, 81, 1 >, Eigen::Matrix< double, 81, 81 > >(const assembler::NonLinearAssemblerData &)> &fun81, const std::function< DScalar2< double, Eigen::Matrix< double, Eigen::Dynamic, 1, 0, SMALL_N, 1 >, Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic, 0, SMALL_N, SMALL_N > >(const assembler::NonLinearAssemblerData &)> &funN, const std::function< DScalar2< double, Eigen::VectorXd, Eigen::MatrixXd >(const assembler::NonLinearAssemblerData &)> &funn)
Definition ElasticityUtils.cpp:180
polyfem::json
nlohmann::json json
Definition Common.hpp:9
polyfem::get_local_disp
void get_local_disp(const assembler::NonLinearAssemblerData &data, const int size, AutoDiffVect &local_disp)
Definition ElasticityUtils.hpp:71
polyfem::gradient_from_energy
Eigen::VectorXd gradient_from_energy(const int size, const int n_bases, const assembler::NonLinearAssemblerData &data, const std::function< DScalar1< double, Eigen::Matrix< double, 6, 1 > >(const assembler::NonLinearAssemblerData &)> &fun6, const std::function< DScalar1< double, Eigen::Matrix< double, 8, 1 > >(const assembler::NonLinearAssemblerData &)> &fun8, const std::function< DScalar1< double, Eigen::Matrix< double, 12, 1 > >(const assembler::NonLinearAssemblerData &)> &fun12, const std::function< DScalar1< double, Eigen::Matrix< double, 18, 1 > >(const assembler::NonLinearAssemblerData &)> &fun18, const std::function< DScalar1< double, Eigen::Matrix< double, 24, 1 > >(const assembler::NonLinearAssemblerData &)> &fun24, const std::function< DScalar1< double, Eigen::Matrix< double, 30, 1 > >(const assembler::NonLinearAssemblerData &)> &fun30, const std::function< DScalar1< double, Eigen::Matrix< double, 60, 1 > >(const assembler::NonLinearAssemblerData &)> &fun60, const std::function< DScalar1< double, Eigen::Matrix< double, 81, 1 > >(const assembler::NonLinearAssemblerData &)> &fun81, const std::function< DScalar1< double, Eigen::Matrix< double, Eigen::Dynamic, 1, 0, SMALL_N, 1 > >(const assembler::NonLinearAssemblerData &)> &funN, const std::function< DScalar1< double, Eigen::Matrix< double, Eigen::Dynamic, 1, 0, BIG_N, 1 > >(const assembler::NonLinearAssemblerData &)> &funBigN, const std::function< DScalar1< double, Eigen::VectorXd >(const assembler::NonLinearAssemblerData &)> &funn)
Definition ElasticityUtils.cpp:81
DScalar2
Automatic differentiation scalar with first- and second-order derivatives.
Definition autodiff.h:493
DiffScalarBase::setVariableCount
static void setVariableCount(size_t value)
Set the independent variable count used by the automatic differentiation layer.
Definition autodiff.h:54