PhysicsBasedAnimationToolkit 0.0.10
Cross-platform C++20 library of algorithms and data structures commonly used in computer graphics research on physically-based simulation.
Loading...
Searching...
No Matches
StableNeoHookeanEnergy.h
Go to the documentation of this file.
1
9
10#ifndef PBAT_PHYSICS_STABLENEOHOOKEANENERGY_H
11#define PBAT_PHYSICS_STABLENEOHOOKEANENERGY_H
12
13#include "pbat/Aliases.h"
14#include "pbat/HostDevice.h"
15#include "pbat/math/linalg/mini/Matrix.h"
16
17#include <cmath>
18
19namespace pbat {
20namespace physics {
21
22template <int Dims>
23struct StableNeoHookeanEnergy;
24
30template <>
31struct StableNeoHookeanEnergy<1>
32{
33 public:
34 template <class TScalar, int M, int N>
35 using SMatrix = pbat::math::linalg::mini::SMatrix<TScalar, M, N>;
36
37 template <class TScalar, int M>
38 using SVector = pbat::math::linalg::mini::SVector<TScalar, M>;
39
40 static auto constexpr kDims = 1;
41
51 template <math::linalg::mini::CReadableVectorizedMatrix TMatrix>
52 PBAT_HOST_DEVICE typename TMatrix::ScalarType
53 eval(TMatrix const& F, typename TMatrix::ScalarType mu, typename TMatrix::ScalarType lambda)
54 const;
55
65 template <math::linalg::mini::CReadableVectorizedMatrix TMatrix>
66 PBAT_HOST_DEVICE SVector<typename TMatrix::ScalarType, 1>
67 grad(TMatrix const& F, typename TMatrix::ScalarType mu, typename TMatrix::ScalarType lambda)
68 const;
69
79 template <math::linalg::mini::CReadableVectorizedMatrix TMatrix>
80 PBAT_HOST_DEVICE SMatrix<typename TMatrix::ScalarType, 1, 1>
81 hessian(TMatrix const& F, typename TMatrix::ScalarType mu, typename TMatrix::ScalarType lambda)
82 const;
83
94 template <
95 math::linalg::mini::CReadableVectorizedMatrix TMatrix,
96 math::linalg::mini::CWriteableVectorizedMatrix TMatrixGF>
97 PBAT_HOST_DEVICE typename TMatrix::ScalarType evalWithGrad(
98 TMatrix const& F,
99 typename TMatrix::ScalarType mu,
100 typename TMatrix::ScalarType lambda,
101 TMatrixGF& gF) const;
102
114 template <
115 math::linalg::mini::CReadableVectorizedMatrix TMatrix,
116 math::linalg::mini::CWriteableVectorizedMatrix TMatrixGF,
117 math::linalg::mini::CWriteableVectorizedMatrix TMatrixHF>
118 PBAT_HOST_DEVICE typename TMatrix::ScalarType evalWithGradAndHessian(
119 TMatrix const& F,
120 typename TMatrix::ScalarType mu,
121 typename TMatrix::ScalarType lambda,
122 TMatrixGF& gF,
123 TMatrixHF& HF) const;
124
135 template <
136 math::linalg::mini::CReadableVectorizedMatrix TMatrix,
137 math::linalg::mini::CWriteableVectorizedMatrix TMatrixGF,
138 math::linalg::mini::CWriteableVectorizedMatrix TMatrixHF>
139 PBAT_HOST_DEVICE void gradAndHessian(
140 TMatrix const& F,
141 typename TMatrix::ScalarType mu,
142 typename TMatrix::ScalarType lambda,
143 TMatrixGF& gF,
144 TMatrixHF& HF) const;
145};
146
147template <math::linalg::mini::CReadableVectorizedMatrix TMatrix>
148PBAT_HOST_DEVICE typename TMatrix::ScalarType StableNeoHookeanEnergy<1>::eval(
149 [[maybe_unused]] TMatrix const& F,
150 [[maybe_unused]] typename TMatrix::ScalarType mu,
151 [[maybe_unused]] typename TMatrix::ScalarType lambda) const
152{
153 using ScalarType = typename TMatrix::ScalarType;
154 ScalarType psi;
155 psi = (1.0 / 2.0) * lambda * ((F[0] - 1 - mu / lambda) * (F[0] - 1 - mu / lambda)) +
156 (1.0 / 2.0) * mu * (((F[0]) * (F[0])) - 1);
157 return psi;
158}
159
169template <math::linalg::mini::CReadableVectorizedMatrix TMatrix>
170PBAT_HOST_DEVICE StableNeoHookeanEnergy<1>::SVector<typename TMatrix::ScalarType, 1>
171StableNeoHookeanEnergy<1>::grad(
172 [[maybe_unused]] TMatrix const& F,
173 [[maybe_unused]] typename TMatrix::ScalarType mu,
174 [[maybe_unused]] typename TMatrix::ScalarType lambda) const
175{
176 using ScalarType = typename TMatrix::ScalarType;
177 SVector<ScalarType, 1> G;
178 G[0] = (1.0 / 2.0) * lambda * (2 * F[0] - 2 - 2 * mu / lambda) + mu * F[0];
179 return G;
180}
181
191template <math::linalg::mini::CReadableVectorizedMatrix TMatrix>
192PBAT_HOST_DEVICE StableNeoHookeanEnergy<1>::SMatrix<typename TMatrix::ScalarType, 1, 1>
193StableNeoHookeanEnergy<1>::hessian(
194 [[maybe_unused]] TMatrix const& F,
195 [[maybe_unused]] typename TMatrix::ScalarType mu,
196 [[maybe_unused]] typename TMatrix::ScalarType lambda) const
197{
198 using ScalarType = typename TMatrix::ScalarType;
199 SMatrix<ScalarType, 1, 1> H;
200 H[0] = lambda + mu;
201 return H;
202}
203
204template <
205 math::linalg::mini::CReadableVectorizedMatrix TMatrix,
206 math::linalg::mini::CWriteableVectorizedMatrix TMatrixGF>
207PBAT_HOST_DEVICE typename TMatrix::ScalarType StableNeoHookeanEnergy<1>::evalWithGrad(
208 [[maybe_unused]] TMatrix const& F,
209 [[maybe_unused]] typename TMatrix::ScalarType mu,
210 [[maybe_unused]] typename TMatrix::ScalarType lambda,
211 TMatrixGF& gF) const
212{
213 static_assert(
214 TMatrixGF::kRows == 1 and TMatrixGF::kCols == 1,
215 "Grad w.r.t. F must have dimensions 1x1");
216 using ScalarType = typename TMatrix::ScalarType;
217 ScalarType psi;
218 ScalarType const a0 = mu / lambda;
219 ScalarType const a1 = (1.0 / 2.0) * lambda;
220 psi = a1 * ((-a0 + F[0] - 1) * (-a0 + F[0] - 1)) + (1.0 / 2.0) * mu * (((F[0]) * (F[0])) - 1);
221 gF[0] = a1 * (-2 * a0 + 2 * F[0] - 2) + mu * F[0];
222 return psi;
223}
224
225template <
226 math::linalg::mini::CReadableVectorizedMatrix TMatrix,
227 math::linalg::mini::CWriteableVectorizedMatrix TMatrixGF,
228 math::linalg::mini::CWriteableVectorizedMatrix TMatrixHF>
229PBAT_HOST_DEVICE typename TMatrix::ScalarType StableNeoHookeanEnergy<1>::evalWithGradAndHessian(
230 [[maybe_unused]] TMatrix const& F,
231 [[maybe_unused]] typename TMatrix::ScalarType mu,
232 [[maybe_unused]] typename TMatrix::ScalarType lambda,
233 TMatrixGF& gF,
234 TMatrixHF& HF) const
235{
236 static_assert(
237 TMatrixGF::kRows == 1 and TMatrixGF::kCols == 1,
238 "Grad w.r.t. F must have dimensions 1x1");
239 static_assert(
240 TMatrixHF::kRows == 1 and TMatrixHF::kCols == 1,
241 "Hessian w.r.t. F must have dimensions 1x1");
242 using ScalarType = typename TMatrix::ScalarType;
243 ScalarType psi;
244 ScalarType const a0 = mu / lambda;
245 ScalarType const a1 = (1.0 / 2.0) * lambda;
246 psi = a1 * ((-a0 + F[0] - 1) * (-a0 + F[0] - 1)) + (1.0 / 2.0) * mu * (((F[0]) * (F[0])) - 1);
247 gF[0] = a1 * (-2 * a0 + 2 * F[0] - 2) + mu * F[0];
248 HF[0] = lambda + mu;
249 return psi;
250}
251
252template <
253 math::linalg::mini::CReadableVectorizedMatrix TMatrix,
254 math::linalg::mini::CWriteableVectorizedMatrix TMatrixGF,
255 math::linalg::mini::CWriteableVectorizedMatrix TMatrixHF>
256PBAT_HOST_DEVICE void StableNeoHookeanEnergy<1>::gradAndHessian(
257 [[maybe_unused]] TMatrix const& F,
258 [[maybe_unused]] typename TMatrix::ScalarType mu,
259 [[maybe_unused]] typename TMatrix::ScalarType lambda,
260 TMatrixGF& gF,
261 TMatrixHF& HF) const
262{
263 static_assert(
264 TMatrixGF::kRows == 1 and TMatrixGF::kCols == 1,
265 "Grad w.r.t. F must have dimensions 1x1");
266 static_assert(
267 TMatrixHF::kRows == 1 and TMatrixHF::kCols == 1,
268 "Hessian w.r.t. F must have dimensions 1x1");
269 gF[0] = (1.0 / 2.0) * lambda * (2 * F[0] - 2 - 2 * mu / lambda) + mu * F[0];
270 HF[0] = lambda + mu;
271}
272
278template <>
279struct StableNeoHookeanEnergy<2>
280{
281 public:
282 template <class TScalar, int M, int N>
283 using SMatrix = pbat::math::linalg::mini::SMatrix<TScalar, M, N>;
284
285 template <class TScalar, int M>
286 using SVector = pbat::math::linalg::mini::SVector<TScalar, M>;
287
288 static auto constexpr kDims = 2;
289
299 template <math::linalg::mini::CReadableVectorizedMatrix TMatrix>
300 PBAT_HOST_DEVICE typename TMatrix::ScalarType
301 eval(TMatrix const& F, typename TMatrix::ScalarType mu, typename TMatrix::ScalarType lambda)
302 const;
303
313 template <math::linalg::mini::CReadableVectorizedMatrix TMatrix>
314 PBAT_HOST_DEVICE SVector<typename TMatrix::ScalarType, 4>
315 grad(TMatrix const& F, typename TMatrix::ScalarType mu, typename TMatrix::ScalarType lambda)
316 const;
317
327 template <math::linalg::mini::CReadableVectorizedMatrix TMatrix>
328 PBAT_HOST_DEVICE SMatrix<typename TMatrix::ScalarType, 4, 4>
329 hessian(TMatrix const& F, typename TMatrix::ScalarType mu, typename TMatrix::ScalarType lambda)
330 const;
331
342 template <
343 math::linalg::mini::CReadableVectorizedMatrix TMatrix,
344 math::linalg::mini::CWriteableVectorizedMatrix TMatrixGF>
345 PBAT_HOST_DEVICE typename TMatrix::ScalarType evalWithGrad(
346 TMatrix const& F,
347 typename TMatrix::ScalarType mu,
348 typename TMatrix::ScalarType lambda,
349 TMatrixGF& gF) const;
350
362 template <
363 math::linalg::mini::CReadableVectorizedMatrix TMatrix,
364 math::linalg::mini::CWriteableVectorizedMatrix TMatrixGF,
365 math::linalg::mini::CWriteableVectorizedMatrix TMatrixHF>
366 PBAT_HOST_DEVICE typename TMatrix::ScalarType evalWithGradAndHessian(
367 TMatrix const& F,
368 typename TMatrix::ScalarType mu,
369 typename TMatrix::ScalarType lambda,
370 TMatrixGF& gF,
371 TMatrixHF& HF) const;
372
383 template <
384 math::linalg::mini::CReadableVectorizedMatrix TMatrix,
385 math::linalg::mini::CWriteableVectorizedMatrix TMatrixGF,
386 math::linalg::mini::CWriteableVectorizedMatrix TMatrixHF>
387 PBAT_HOST_DEVICE void gradAndHessian(
388 TMatrix const& F,
389 typename TMatrix::ScalarType mu,
390 typename TMatrix::ScalarType lambda,
391 TMatrixGF& gF,
392 TMatrixHF& HF) const;
393};
394
395template <math::linalg::mini::CReadableVectorizedMatrix TMatrix>
396PBAT_HOST_DEVICE typename TMatrix::ScalarType StableNeoHookeanEnergy<2>::eval(
397 [[maybe_unused]] TMatrix const& F,
398 [[maybe_unused]] typename TMatrix::ScalarType mu,
399 [[maybe_unused]] typename TMatrix::ScalarType lambda) const
400{
401 using ScalarType = typename TMatrix::ScalarType;
402 ScalarType psi;
403 psi = (1.0 / 2.0) * lambda *
404 ((F[0] * F[3] - F[1] * F[2] - 1 - mu / lambda) *
405 (F[0] * F[3] - F[1] * F[2] - 1 - mu / lambda)) +
406 (1.0 / 2.0) * mu *
407 (((F[0]) * (F[0])) + ((F[1]) * (F[1])) + ((F[2]) * (F[2])) + ((F[3]) * (F[3])) - 2);
408 return psi;
409}
410
420template <math::linalg::mini::CReadableVectorizedMatrix TMatrix>
421PBAT_HOST_DEVICE StableNeoHookeanEnergy<2>::SVector<typename TMatrix::ScalarType, 4>
422StableNeoHookeanEnergy<2>::grad(
423 [[maybe_unused]] TMatrix const& F,
424 [[maybe_unused]] typename TMatrix::ScalarType mu,
425 [[maybe_unused]] typename TMatrix::ScalarType lambda) const
426{
427 using ScalarType = typename TMatrix::ScalarType;
428 SVector<ScalarType, 4> G;
429 ScalarType const a0 = lambda * (F[0] * F[3] - F[1] * F[2] - 1 - mu / lambda);
430 G[0] = a0 * F[3] + mu * F[0];
431 G[1] = -a0 * F[2] + mu * F[1];
432 G[2] = -a0 * F[1] + mu * F[2];
433 G[3] = a0 * F[0] + mu * F[3];
434 return G;
435}
436
446template <math::linalg::mini::CReadableVectorizedMatrix TMatrix>
447PBAT_HOST_DEVICE StableNeoHookeanEnergy<2>::SMatrix<typename TMatrix::ScalarType, 4, 4>
448StableNeoHookeanEnergy<2>::hessian(
449 [[maybe_unused]] TMatrix const& F,
450 [[maybe_unused]] typename TMatrix::ScalarType mu,
451 [[maybe_unused]] typename TMatrix::ScalarType lambda) const
452{
453 using ScalarType = typename TMatrix::ScalarType;
454 SMatrix<ScalarType, 4, 4> H;
455 ScalarType const a0 = lambda * F[3];
456 ScalarType const a1 = -a0 * F[2];
457 ScalarType const a2 = -a0 * F[1];
458 ScalarType const a3 = lambda * (F[0] * F[3] - F[1] * F[2] - 1 - mu / lambda);
459 ScalarType const a4 = a3 + lambda * F[0] * F[3];
460 ScalarType const a5 = -a3 + lambda * F[1] * F[2];
461 ScalarType const a6 = lambda * F[0];
462 ScalarType const a7 = -a6 * F[2];
463 ScalarType const a8 = -a6 * F[1];
464 H[0] = lambda * ((F[3]) * (F[3])) + mu;
465 H[1] = a1;
466 H[2] = a2;
467 H[3] = a4;
468 H[4] = a1;
469 H[5] = lambda * ((F[2]) * (F[2])) + mu;
470 H[6] = a5;
471 H[7] = a7;
472 H[8] = a2;
473 H[9] = a5;
474 H[10] = lambda * ((F[1]) * (F[1])) + mu;
475 H[11] = a8;
476 H[12] = a4;
477 H[13] = a7;
478 H[14] = a8;
479 H[15] = lambda * ((F[0]) * (F[0])) + mu;
480 return H;
481}
482
483template <
484 math::linalg::mini::CReadableVectorizedMatrix TMatrix,
485 math::linalg::mini::CWriteableVectorizedMatrix TMatrixGF>
486PBAT_HOST_DEVICE typename TMatrix::ScalarType StableNeoHookeanEnergy<2>::evalWithGrad(
487 [[maybe_unused]] TMatrix const& F,
488 [[maybe_unused]] typename TMatrix::ScalarType mu,
489 [[maybe_unused]] typename TMatrix::ScalarType lambda,
490 TMatrixGF& gF) const
491{
492 static_assert(
493 TMatrixGF::kRows == 4 and TMatrixGF::kCols == 1,
494 "Grad w.r.t. F must have dimensions 4x1");
495 using ScalarType = typename TMatrix::ScalarType;
496 ScalarType psi;
497 ScalarType const a0 = F[0] * F[3] - F[1] * F[2] - 1 - mu / lambda;
498 ScalarType const a1 = a0 * lambda;
499 psi = (1.0 / 2.0) * ((a0) * (a0)) * lambda +
500 (1.0 / 2.0) * mu *
501 (((F[0]) * (F[0])) + ((F[1]) * (F[1])) + ((F[2]) * (F[2])) + ((F[3]) * (F[3])) - 2);
502 gF[0] = a1 * F[3] + mu * F[0];
503 gF[1] = -a1 * F[2] + mu * F[1];
504 gF[2] = -a1 * F[1] + mu * F[2];
505 gF[3] = a1 * F[0] + mu * F[3];
506 return psi;
507}
508
509template <
510 math::linalg::mini::CReadableVectorizedMatrix TMatrix,
511 math::linalg::mini::CWriteableVectorizedMatrix TMatrixGF,
512 math::linalg::mini::CWriteableVectorizedMatrix TMatrixHF>
513PBAT_HOST_DEVICE typename TMatrix::ScalarType StableNeoHookeanEnergy<2>::evalWithGradAndHessian(
514 [[maybe_unused]] TMatrix const& F,
515 [[maybe_unused]] typename TMatrix::ScalarType mu,
516 [[maybe_unused]] typename TMatrix::ScalarType lambda,
517 TMatrixGF& gF,
518 TMatrixHF& HF) const
519{
520 static_assert(
521 TMatrixGF::kRows == 4 and TMatrixGF::kCols == 1,
522 "Grad w.r.t. F must have dimensions 4x1");
523 static_assert(
524 TMatrixHF::kRows == 4 and TMatrixHF::kCols == 4,
525 "Hessian w.r.t. F must have dimensions 4x4");
526 using ScalarType = typename TMatrix::ScalarType;
527 ScalarType psi;
528 ScalarType const a0 = ((F[0]) * (F[0]));
529 ScalarType const a1 = ((F[1]) * (F[1]));
530 ScalarType const a2 = ((F[2]) * (F[2]));
531 ScalarType const a3 = ((F[3]) * (F[3]));
532 ScalarType const a4 = F[0] * F[3] - F[1] * F[2] - 1 - mu / lambda;
533 ScalarType const a5 = a4 * lambda;
534 ScalarType const a6 = lambda * F[3];
535 ScalarType const a7 = -a6 * F[2];
536 ScalarType const a8 = -a6 * F[1];
537 ScalarType const a9 = a5 + lambda * F[0] * F[3];
538 ScalarType const a10 = -a5 + lambda * F[1] * F[2];
539 ScalarType const a11 = lambda * F[0];
540 ScalarType const a12 = -a11 * F[2];
541 ScalarType const a13 = -a11 * F[1];
542 psi = (1.0 / 2.0) * ((a4) * (a4)) * lambda + (1.0 / 2.0) * mu * (a0 + a1 + a2 + a3 - 2);
543 gF[0] = a5 * F[3] + mu * F[0];
544 gF[1] = -a5 * F[2] + mu * F[1];
545 gF[2] = -a5 * F[1] + mu * F[2];
546 gF[3] = a5 * F[0] + mu * F[3];
547 HF[0] = a3 * lambda + mu;
548 HF[1] = a7;
549 HF[2] = a8;
550 HF[3] = a9;
551 HF[4] = a7;
552 HF[5] = a2 * lambda + mu;
553 HF[6] = a10;
554 HF[7] = a12;
555 HF[8] = a8;
556 HF[9] = a10;
557 HF[10] = a1 * lambda + mu;
558 HF[11] = a13;
559 HF[12] = a9;
560 HF[13] = a12;
561 HF[14] = a13;
562 HF[15] = a0 * lambda + mu;
563 return psi;
564}
565
566template <
567 math::linalg::mini::CReadableVectorizedMatrix TMatrix,
568 math::linalg::mini::CWriteableVectorizedMatrix TMatrixGF,
569 math::linalg::mini::CWriteableVectorizedMatrix TMatrixHF>
570PBAT_HOST_DEVICE void StableNeoHookeanEnergy<2>::gradAndHessian(
571 [[maybe_unused]] TMatrix const& F,
572 [[maybe_unused]] typename TMatrix::ScalarType mu,
573 [[maybe_unused]] typename TMatrix::ScalarType lambda,
574 TMatrixGF& gF,
575 TMatrixHF& HF) const
576{
577 static_assert(
578 TMatrixGF::kRows == 4 and TMatrixGF::kCols == 1,
579 "Grad w.r.t. F must have dimensions 4x1");
580 static_assert(
581 TMatrixHF::kRows == 4 and TMatrixHF::kCols == 4,
582 "Hessian w.r.t. F must have dimensions 4x4");
583 using ScalarType = typename TMatrix::ScalarType;
584 ScalarType const a0 = lambda * (F[0] * F[3] - F[1] * F[2] - 1 - mu / lambda);
585 ScalarType const a1 = lambda * F[3];
586 ScalarType const a2 = -a1 * F[2];
587 ScalarType const a3 = -a1 * F[1];
588 ScalarType const a4 = a0 + lambda * F[0] * F[3];
589 ScalarType const a5 = -a0 + lambda * F[1] * F[2];
590 ScalarType const a6 = lambda * F[0];
591 ScalarType const a7 = -a6 * F[2];
592 ScalarType const a8 = -a6 * F[1];
593 gF[0] = a0 * F[3] + mu * F[0];
594 gF[1] = -a0 * F[2] + mu * F[1];
595 gF[2] = -a0 * F[1] + mu * F[2];
596 gF[3] = a0 * F[0] + mu * F[3];
597 HF[0] = lambda * ((F[3]) * (F[3])) + mu;
598 HF[1] = a2;
599 HF[2] = a3;
600 HF[3] = a4;
601 HF[4] = a2;
602 HF[5] = lambda * ((F[2]) * (F[2])) + mu;
603 HF[6] = a5;
604 HF[7] = a7;
605 HF[8] = a3;
606 HF[9] = a5;
607 HF[10] = lambda * ((F[1]) * (F[1])) + mu;
608 HF[11] = a8;
609 HF[12] = a4;
610 HF[13] = a7;
611 HF[14] = a8;
612 HF[15] = lambda * ((F[0]) * (F[0])) + mu;
613}
614
620template <>
621struct StableNeoHookeanEnergy<3>
622{
623 public:
624 template <class TScalar, int M, int N>
625 using SMatrix = pbat::math::linalg::mini::SMatrix<TScalar, M, N>;
626
627 template <class TScalar, int M>
628 using SVector = pbat::math::linalg::mini::SVector<TScalar, M>;
629
630 static auto constexpr kDims = 3;
631
641 template <math::linalg::mini::CReadableVectorizedMatrix TMatrix>
642 PBAT_HOST_DEVICE typename TMatrix::ScalarType
643 eval(TMatrix const& F, typename TMatrix::ScalarType mu, typename TMatrix::ScalarType lambda)
644 const;
645
655 template <math::linalg::mini::CReadableVectorizedMatrix TMatrix>
656 PBAT_HOST_DEVICE SVector<typename TMatrix::ScalarType, 9>
657 grad(TMatrix const& F, typename TMatrix::ScalarType mu, typename TMatrix::ScalarType lambda)
658 const;
659
669 template <math::linalg::mini::CReadableVectorizedMatrix TMatrix>
670 PBAT_HOST_DEVICE SMatrix<typename TMatrix::ScalarType, 9, 9>
671 hessian(TMatrix const& F, typename TMatrix::ScalarType mu, typename TMatrix::ScalarType lambda)
672 const;
673
684 template <
685 math::linalg::mini::CReadableVectorizedMatrix TMatrix,
686 math::linalg::mini::CWriteableVectorizedMatrix TMatrixGF>
687 PBAT_HOST_DEVICE typename TMatrix::ScalarType evalWithGrad(
688 TMatrix const& F,
689 typename TMatrix::ScalarType mu,
690 typename TMatrix::ScalarType lambda,
691 TMatrixGF& gF) const;
692
704 template <
705 math::linalg::mini::CReadableVectorizedMatrix TMatrix,
706 math::linalg::mini::CWriteableVectorizedMatrix TMatrixGF,
707 math::linalg::mini::CWriteableVectorizedMatrix TMatrixHF>
708 PBAT_HOST_DEVICE typename TMatrix::ScalarType evalWithGradAndHessian(
709 TMatrix const& F,
710 typename TMatrix::ScalarType mu,
711 typename TMatrix::ScalarType lambda,
712 TMatrixGF& gF,
713 TMatrixHF& HF) const;
714
725 template <
726 math::linalg::mini::CReadableVectorizedMatrix TMatrix,
727 math::linalg::mini::CWriteableVectorizedMatrix TMatrixGF,
728 math::linalg::mini::CWriteableVectorizedMatrix TMatrixHF>
729 PBAT_HOST_DEVICE void gradAndHessian(
730 TMatrix const& F,
731 typename TMatrix::ScalarType mu,
732 typename TMatrix::ScalarType lambda,
733 TMatrixGF& gF,
734 TMatrixHF& HF) const;
735};
736
737template <math::linalg::mini::CReadableVectorizedMatrix TMatrix>
738PBAT_HOST_DEVICE typename TMatrix::ScalarType StableNeoHookeanEnergy<3>::eval(
739 [[maybe_unused]] TMatrix const& F,
740 [[maybe_unused]] typename TMatrix::ScalarType mu,
741 [[maybe_unused]] typename TMatrix::ScalarType lambda) const
742{
743 using ScalarType = typename TMatrix::ScalarType;
744 ScalarType psi;
745 psi = (1.0 / 2.0) * lambda *
746 ((F[0] * F[4] * F[8] - F[0] * F[5] * F[7] - F[1] * F[3] * F[8] + F[1] * F[5] * F[6] +
747 F[2] * F[3] * F[7] - F[2] * F[4] * F[6] - 1 - mu / lambda) *
748 (F[0] * F[4] * F[8] - F[0] * F[5] * F[7] - F[1] * F[3] * F[8] + F[1] * F[5] * F[6] +
749 F[2] * F[3] * F[7] - F[2] * F[4] * F[6] - 1 - mu / lambda)) +
750 (1.0 / 2.0) * mu *
751 (((F[0]) * (F[0])) + ((F[1]) * (F[1])) + ((F[2]) * (F[2])) + ((F[3]) * (F[3])) +
752 ((F[4]) * (F[4])) + ((F[5]) * (F[5])) + ((F[6]) * (F[6])) + ((F[7]) * (F[7])) +
753 ((F[8]) * (F[8])) - 3);
754 return psi;
755}
756
766template <math::linalg::mini::CReadableVectorizedMatrix TMatrix>
767PBAT_HOST_DEVICE StableNeoHookeanEnergy<3>::SVector<typename TMatrix::ScalarType, 9>
768StableNeoHookeanEnergy<3>::grad(
769 [[maybe_unused]] TMatrix const& F,
770 [[maybe_unused]] typename TMatrix::ScalarType mu,
771 [[maybe_unused]] typename TMatrix::ScalarType lambda) const
772{
773 using ScalarType = typename TMatrix::ScalarType;
774 SVector<ScalarType, 9> G;
775 ScalarType const a0 = F[5] * F[7];
776 ScalarType const a1 = F[3] * F[8];
777 ScalarType const a2 = F[4] * F[6];
778 ScalarType const a3 = (1.0 / 2.0) * lambda *
779 (-a0 * F[0] - a1 * F[1] - a2 * F[2] + F[0] * F[4] * F[8] +
780 F[1] * F[5] * F[6] + F[2] * F[3] * F[7] - 1 - mu / lambda);
781 ScalarType const a4 = 2 * F[8];
782 ScalarType const a5 = 2 * F[2];
783 ScalarType const a6 = 2 * F[0];
784 ScalarType const a7 = 2 * F[1];
785 G[0] = a3 * (-2 * a0 + 2 * F[4] * F[8]) + mu * F[0];
786 G[1] = a3 * (-2 * a1 + 2 * F[5] * F[6]) + mu * F[1];
787 G[2] = a3 * (-2 * a2 + 2 * F[3] * F[7]) + mu * F[2];
788 G[3] = a3 * (-a4 * F[1] + 2 * F[2] * F[7]) + mu * F[3];
789 G[4] = a3 * (a4 * F[0] - a5 * F[6]) + mu * F[4];
790 G[5] = a3 * (-a6 * F[7] + 2 * F[1] * F[6]) + mu * F[5];
791 G[6] = a3 * (-a5 * F[4] + a7 * F[5]) + mu * F[6];
792 G[7] = a3 * (-a6 * F[5] + 2 * F[2] * F[3]) + mu * F[7];
793 G[8] = a3 * (a6 * F[4] - a7 * F[3]) + mu * F[8];
794 return G;
795}
796
806template <math::linalg::mini::CReadableVectorizedMatrix TMatrix>
807PBAT_HOST_DEVICE StableNeoHookeanEnergy<3>::SMatrix<typename TMatrix::ScalarType, 9, 9>
808StableNeoHookeanEnergy<3>::hessian(
809 [[maybe_unused]] TMatrix const& F,
810 [[maybe_unused]] typename TMatrix::ScalarType mu,
811 [[maybe_unused]] typename TMatrix::ScalarType lambda) const
812{
813 using ScalarType = typename TMatrix::ScalarType;
814 SMatrix<ScalarType, 9, 9> H;
815 ScalarType const a0 = F[4] * F[8];
816 ScalarType const a1 = F[5] * F[7];
817 ScalarType const a2 = a0 - a1;
818 ScalarType const a3 = (1.0 / 2.0) * lambda;
819 ScalarType const a4 = a3 * (2 * a0 - 2 * a1);
820 ScalarType const a5 = F[3] * F[8];
821 ScalarType const a6 = -a5 + F[5] * F[6];
822 ScalarType const a7 = F[3] * F[7];
823 ScalarType const a8 = F[4] * F[6];
824 ScalarType const a9 = a7 - a8;
825 ScalarType const a10 = F[1] * F[8];
826 ScalarType const a11 = -a10 + F[2] * F[7];
827 ScalarType const a12 = F[0] * F[8];
828 ScalarType const a13 = F[2] * F[6];
829 ScalarType const a14 = a12 - a13;
830 ScalarType const a15 = lambda * (-a1 * F[0] - a5 * F[1] - a8 * F[2] + F[0] * F[4] * F[8] +
831 F[1] * F[5] * F[6] + F[2] * F[3] * F[7] - 1 - mu / lambda);
832 ScalarType const a16 = a15 * F[8];
833 ScalarType const a17 = F[0] * F[7];
834 ScalarType const a18 = -a17 + F[1] * F[6];
835 ScalarType const a19 = a15 * F[7];
836 ScalarType const a20 = -a19;
837 ScalarType const a21 = F[1] * F[5];
838 ScalarType const a22 = F[2] * F[4];
839 ScalarType const a23 = a21 - a22;
840 ScalarType const a24 = F[0] * F[5];
841 ScalarType const a25 = -a24 + F[2] * F[3];
842 ScalarType const a26 = a15 * F[5];
843 ScalarType const a27 = -a26;
844 ScalarType const a28 = F[0] * F[4];
845 ScalarType const a29 = F[1] * F[3];
846 ScalarType const a30 = a28 - a29;
847 ScalarType const a31 = a15 * F[4];
848 ScalarType const a32 = a3 * (-2 * a5 + 2 * F[5] * F[6]);
849 ScalarType const a33 = -a16;
850 ScalarType const a34 = a15 * F[6];
851 ScalarType const a35 = a15 * F[3];
852 ScalarType const a36 = -a35;
853 ScalarType const a37 = a3 * (2 * a7 - 2 * a8);
854 ScalarType const a38 = -a34;
855 ScalarType const a39 = -a31;
856 ScalarType const a40 = a3 * (-2 * a10 + 2 * F[2] * F[7]);
857 ScalarType const a41 = a15 * F[2];
858 ScalarType const a42 = a15 * F[1];
859 ScalarType const a43 = -a42;
860 ScalarType const a44 = a3 * (2 * a12 - 2 * a13);
861 ScalarType const a45 = -a41;
862 ScalarType const a46 = a15 * F[0];
863 ScalarType const a47 = a3 * (-2 * a17 + 2 * F[1] * F[6]);
864 ScalarType const a48 = -a46;
865 ScalarType const a49 = a3 * (2 * a21 - 2 * a22);
866 ScalarType const a50 = a3 * (-2 * a24 + 2 * F[2] * F[3]);
867 ScalarType const a51 = a3 * (2 * a28 - 2 * a29);
868 H[0] = a2 * a4 + mu;
869 H[1] = a4 * a6;
870 H[2] = a4 * a9;
871 H[3] = a11 * a4;
872 H[4] = a14 * a4 + a16;
873 H[5] = a18 * a4 + a20;
874 H[6] = a23 * a4;
875 H[7] = a25 * a4 + a27;
876 H[8] = a30 * a4 + a31;
877 H[9] = a2 * a32;
878 H[10] = a32 * a6 + mu;
879 H[11] = a32 * a9;
880 H[12] = a11 * a32 + a33;
881 H[13] = a14 * a32;
882 H[14] = a18 * a32 + a34;
883 H[15] = a23 * a32 + a26;
884 H[16] = a25 * a32;
885 H[17] = a30 * a32 + a36;
886 H[18] = a2 * a37;
887 H[19] = a37 * a6;
888 H[20] = a37 * a9 + mu;
889 H[21] = a11 * a37 + a19;
890 H[22] = a14 * a37 + a38;
891 H[23] = a18 * a37;
892 H[24] = a23 * a37 + a39;
893 H[25] = a25 * a37 + a35;
894 H[26] = a30 * a37;
895 H[27] = a2 * a40;
896 H[28] = a33 + a40 * a6;
897 H[29] = a19 + a40 * a9;
898 H[30] = a11 * a40 + mu;
899 H[31] = a14 * a40;
900 H[32] = a18 * a40;
901 H[33] = a23 * a40;
902 H[34] = a25 * a40 + a41;
903 H[35] = a30 * a40 + a43;
904 H[36] = a16 + a2 * a44;
905 H[37] = a44 * a6;
906 H[38] = a38 + a44 * a9;
907 H[39] = a11 * a44;
908 H[40] = a14 * a44 + mu;
909 H[41] = a18 * a44;
910 H[42] = a23 * a44 + a45;
911 H[43] = a25 * a44;
912 H[44] = a30 * a44 + a46;
913 H[45] = a2 * a47 + a20;
914 H[46] = a34 + a47 * a6;
915 H[47] = a47 * a9;
916 H[48] = a11 * a47;
917 H[49] = a14 * a47;
918 H[50] = a18 * a47 + mu;
919 H[51] = a23 * a47 + a42;
920 H[52] = a25 * a47 + a48;
921 H[53] = a30 * a47;
922 H[54] = a2 * a49;
923 H[55] = a26 + a49 * a6;
924 H[56] = a39 + a49 * a9;
925 H[57] = a11 * a49;
926 H[58] = a14 * a49 + a45;
927 H[59] = a18 * a49 + a42;
928 H[60] = a23 * a49 + mu;
929 H[61] = a25 * a49;
930 H[62] = a30 * a49;
931 H[63] = a2 * a50 + a27;
932 H[64] = a50 * a6;
933 H[65] = a35 + a50 * a9;
934 H[66] = a11 * a50 + a41;
935 H[67] = a14 * a50;
936 H[68] = a18 * a50 + a48;
937 H[69] = a23 * a50;
938 H[70] = a25 * a50 + mu;
939 H[71] = a30 * a50;
940 H[72] = a2 * a51 + a31;
941 H[73] = a36 + a51 * a6;
942 H[74] = a51 * a9;
943 H[75] = a11 * a51 + a43;
944 H[76] = a14 * a51 + a46;
945 H[77] = a18 * a51;
946 H[78] = a23 * a51;
947 H[79] = a25 * a51;
948 H[80] = a30 * a51 + mu;
949 return H;
950}
951
952template <
953 math::linalg::mini::CReadableVectorizedMatrix TMatrix,
954 math::linalg::mini::CWriteableVectorizedMatrix TMatrixGF>
955PBAT_HOST_DEVICE typename TMatrix::ScalarType StableNeoHookeanEnergy<3>::evalWithGrad(
956 [[maybe_unused]] TMatrix const& F,
957 [[maybe_unused]] typename TMatrix::ScalarType mu,
958 [[maybe_unused]] typename TMatrix::ScalarType lambda,
959 TMatrixGF& gF) const
960{
961 static_assert(
962 TMatrixGF::kRows == 9 and TMatrixGF::kCols == 1,
963 "Grad w.r.t. F must have dimensions 9x1");
964 using ScalarType = typename TMatrix::ScalarType;
965 ScalarType psi;
966 ScalarType const a0 = F[5] * F[7];
967 ScalarType const a1 = F[3] * F[8];
968 ScalarType const a2 = F[4] * F[6];
969 ScalarType const a3 = -a0 * F[0] - a1 * F[1] - a2 * F[2] + F[0] * F[4] * F[8] +
970 F[1] * F[5] * F[6] + F[2] * F[3] * F[7] - 1 - mu / lambda;
971 ScalarType const a4 = (1.0 / 2.0) * lambda;
972 ScalarType const a5 = a3 * a4;
973 ScalarType const a6 = 2 * F[8];
974 ScalarType const a7 = 2 * F[2];
975 ScalarType const a8 = 2 * F[0];
976 ScalarType const a9 = 2 * F[1];
977 psi = ((a3) * (a3)) * a4 + (1.0 / 2.0) * mu *
978 (((F[0]) * (F[0])) + ((F[1]) * (F[1])) + ((F[2]) * (F[2])) +
979 ((F[3]) * (F[3])) + ((F[4]) * (F[4])) + ((F[5]) * (F[5])) +
980 ((F[6]) * (F[6])) + ((F[7]) * (F[7])) + ((F[8]) * (F[8])) - 3);
981 gF[0] = a5 * (-2 * a0 + 2 * F[4] * F[8]) + mu * F[0];
982 gF[1] = a5 * (-2 * a1 + 2 * F[5] * F[6]) + mu * F[1];
983 gF[2] = a5 * (-2 * a2 + 2 * F[3] * F[7]) + mu * F[2];
984 gF[3] = a5 * (-a6 * F[1] + 2 * F[2] * F[7]) + mu * F[3];
985 gF[4] = a5 * (a6 * F[0] - a7 * F[6]) + mu * F[4];
986 gF[5] = a5 * (-a8 * F[7] + 2 * F[1] * F[6]) + mu * F[5];
987 gF[6] = a5 * (-a7 * F[4] + a9 * F[5]) + mu * F[6];
988 gF[7] = a5 * (-a8 * F[5] + 2 * F[2] * F[3]) + mu * F[7];
989 gF[8] = a5 * (a8 * F[4] - a9 * F[3]) + mu * F[8];
990 return psi;
991}
992
993template <
994 math::linalg::mini::CReadableVectorizedMatrix TMatrix,
995 math::linalg::mini::CWriteableVectorizedMatrix TMatrixGF,
996 math::linalg::mini::CWriteableVectorizedMatrix TMatrixHF>
997PBAT_HOST_DEVICE typename TMatrix::ScalarType StableNeoHookeanEnergy<3>::evalWithGradAndHessian(
998 [[maybe_unused]] TMatrix const& F,
999 [[maybe_unused]] typename TMatrix::ScalarType mu,
1000 [[maybe_unused]] typename TMatrix::ScalarType lambda,
1001 TMatrixGF& gF,
1002 TMatrixHF& HF) const
1003{
1004 static_assert(
1005 TMatrixGF::kRows == 9 and TMatrixGF::kCols == 1,
1006 "Grad w.r.t. F must have dimensions 9x1");
1007 static_assert(
1008 TMatrixHF::kRows == 9 and TMatrixHF::kCols == 9,
1009 "Hessian w.r.t. F must have dimensions 9x9");
1010 using ScalarType = typename TMatrix::ScalarType;
1011 ScalarType psi;
1012 ScalarType const a0 = F[5] * F[7];
1013 ScalarType const a1 = F[3] * F[8];
1014 ScalarType const a2 = F[4] * F[6];
1015 ScalarType const a3 = -a0 * F[0] - a1 * F[1] - a2 * F[2] + F[0] * F[4] * F[8] +
1016 F[1] * F[5] * F[6] + F[2] * F[3] * F[7] - 1 - mu / lambda;
1017 ScalarType const a4 = (1.0 / 2.0) * lambda;
1018 ScalarType const a5 = F[4] * F[8];
1019 ScalarType const a6 = -2 * a0 + 2 * a5;
1020 ScalarType const a7 = a3 * a4;
1021 ScalarType const a8 = -2 * a1 + 2 * F[5] * F[6];
1022 ScalarType const a9 = F[3] * F[7];
1023 ScalarType const a10 = -2 * a2 + 2 * a9;
1024 ScalarType const a11 = F[1] * F[8];
1025 ScalarType const a12 = -2 * a11 + 2 * F[2] * F[7];
1026 ScalarType const a13 = F[0] * F[8];
1027 ScalarType const a14 = F[2] * F[6];
1028 ScalarType const a15 = 2 * a13 - 2 * a14;
1029 ScalarType const a16 = F[0] * F[7];
1030 ScalarType const a17 = -2 * a16 + 2 * F[1] * F[6];
1031 ScalarType const a18 = F[1] * F[5];
1032 ScalarType const a19 = F[2] * F[4];
1033 ScalarType const a20 = 2 * a18 - 2 * a19;
1034 ScalarType const a21 = F[0] * F[5];
1035 ScalarType const a22 = -2 * a21 + 2 * F[2] * F[3];
1036 ScalarType const a23 = F[0] * F[4];
1037 ScalarType const a24 = F[1] * F[3];
1038 ScalarType const a25 = 2 * a23 - 2 * a24;
1039 ScalarType const a26 = a4 * (-a0 + a5);
1040 ScalarType const a27 = a3 * lambda;
1041 ScalarType const a28 = a27 * F[8];
1042 ScalarType const a29 = a27 * F[7];
1043 ScalarType const a30 = -a29;
1044 ScalarType const a31 = a27 * F[5];
1045 ScalarType const a32 = -a31;
1046 ScalarType const a33 = a27 * F[4];
1047 ScalarType const a34 = a4 * (-a1 + F[5] * F[6]);
1048 ScalarType const a35 = -a28;
1049 ScalarType const a36 = a27 * F[6];
1050 ScalarType const a37 = a27 * F[3];
1051 ScalarType const a38 = -a37;
1052 ScalarType const a39 = a4 * (-a2 + a9);
1053 ScalarType const a40 = -a36;
1054 ScalarType const a41 = -a33;
1055 ScalarType const a42 = a4 * (-a11 + F[2] * F[7]);
1056 ScalarType const a43 = a27 * F[2];
1057 ScalarType const a44 = a27 * F[1];
1058 ScalarType const a45 = -a44;
1059 ScalarType const a46 = a4 * (a13 - a14);
1060 ScalarType const a47 = -a43;
1061 ScalarType const a48 = a27 * F[0];
1062 ScalarType const a49 = a4 * (-a16 + F[1] * F[6]);
1063 ScalarType const a50 = -a48;
1064 ScalarType const a51 = a4 * (a18 - a19);
1065 ScalarType const a52 = a4 * (-a21 + F[2] * F[3]);
1066 ScalarType const a53 = a4 * (a23 - a24);
1067 psi = ((a3) * (a3)) * a4 + (1.0 / 2.0) * mu *
1068 (((F[0]) * (F[0])) + ((F[1]) * (F[1])) + ((F[2]) * (F[2])) +
1069 ((F[3]) * (F[3])) + ((F[4]) * (F[4])) + ((F[5]) * (F[5])) +
1070 ((F[6]) * (F[6])) + ((F[7]) * (F[7])) + ((F[8]) * (F[8])) - 3);
1071 gF[0] = a6 * a7 + mu * F[0];
1072 gF[1] = a7 * a8 + mu * F[1];
1073 gF[2] = a10 * a7 + mu * F[2];
1074 gF[3] = a12 * a7 + mu * F[3];
1075 gF[4] = a15 * a7 + mu * F[4];
1076 gF[5] = a17 * a7 + mu * F[5];
1077 gF[6] = a20 * a7 + mu * F[6];
1078 gF[7] = a22 * a7 + mu * F[7];
1079 gF[8] = a25 * a7 + mu * F[8];
1080 HF[0] = a26 * a6 + mu;
1081 HF[1] = a26 * a8;
1082 HF[2] = a10 * a26;
1083 HF[3] = a12 * a26;
1084 HF[4] = a15 * a26 + a28;
1085 HF[5] = a17 * a26 + a30;
1086 HF[6] = a20 * a26;
1087 HF[7] = a22 * a26 + a32;
1088 HF[8] = a25 * a26 + a33;
1089 HF[9] = a34 * a6;
1090 HF[10] = a34 * a8 + mu;
1091 HF[11] = a10 * a34;
1092 HF[12] = a12 * a34 + a35;
1093 HF[13] = a15 * a34;
1094 HF[14] = a17 * a34 + a36;
1095 HF[15] = a20 * a34 + a31;
1096 HF[16] = a22 * a34;
1097 HF[17] = a25 * a34 + a38;
1098 HF[18] = a39 * a6;
1099 HF[19] = a39 * a8;
1100 HF[20] = a10 * a39 + mu;
1101 HF[21] = a12 * a39 + a29;
1102 HF[22] = a15 * a39 + a40;
1103 HF[23] = a17 * a39;
1104 HF[24] = a20 * a39 + a41;
1105 HF[25] = a22 * a39 + a37;
1106 HF[26] = a25 * a39;
1107 HF[27] = a42 * a6;
1108 HF[28] = a35 + a42 * a8;
1109 HF[29] = a10 * a42 + a29;
1110 HF[30] = a12 * a42 + mu;
1111 HF[31] = a15 * a42;
1112 HF[32] = a17 * a42;
1113 HF[33] = a20 * a42;
1114 HF[34] = a22 * a42 + a43;
1115 HF[35] = a25 * a42 + a45;
1116 HF[36] = a28 + a46 * a6;
1117 HF[37] = a46 * a8;
1118 HF[38] = a10 * a46 + a40;
1119 HF[39] = a12 * a46;
1120 HF[40] = a15 * a46 + mu;
1121 HF[41] = a17 * a46;
1122 HF[42] = a20 * a46 + a47;
1123 HF[43] = a22 * a46;
1124 HF[44] = a25 * a46 + a48;
1125 HF[45] = a30 + a49 * a6;
1126 HF[46] = a36 + a49 * a8;
1127 HF[47] = a10 * a49;
1128 HF[48] = a12 * a49;
1129 HF[49] = a15 * a49;
1130 HF[50] = a17 * a49 + mu;
1131 HF[51] = a20 * a49 + a44;
1132 HF[52] = a22 * a49 + a50;
1133 HF[53] = a25 * a49;
1134 HF[54] = a51 * a6;
1135 HF[55] = a31 + a51 * a8;
1136 HF[56] = a10 * a51 + a41;
1137 HF[57] = a12 * a51;
1138 HF[58] = a15 * a51 + a47;
1139 HF[59] = a17 * a51 + a44;
1140 HF[60] = a20 * a51 + mu;
1141 HF[61] = a22 * a51;
1142 HF[62] = a25 * a51;
1143 HF[63] = a32 + a52 * a6;
1144 HF[64] = a52 * a8;
1145 HF[65] = a10 * a52 + a37;
1146 HF[66] = a12 * a52 + a43;
1147 HF[67] = a15 * a52;
1148 HF[68] = a17 * a52 + a50;
1149 HF[69] = a20 * a52;
1150 HF[70] = a22 * a52 + mu;
1151 HF[71] = a25 * a52;
1152 HF[72] = a33 + a53 * a6;
1153 HF[73] = a38 + a53 * a8;
1154 HF[74] = a10 * a53;
1155 HF[75] = a12 * a53 + a45;
1156 HF[76] = a15 * a53 + a48;
1157 HF[77] = a17 * a53;
1158 HF[78] = a20 * a53;
1159 HF[79] = a22 * a53;
1160 HF[80] = a25 * a53 + mu;
1161 return psi;
1162}
1163
1164template <
1165 math::linalg::mini::CReadableVectorizedMatrix TMatrix,
1166 math::linalg::mini::CWriteableVectorizedMatrix TMatrixGF,
1167 math::linalg::mini::CWriteableVectorizedMatrix TMatrixHF>
1168PBAT_HOST_DEVICE void StableNeoHookeanEnergy<3>::gradAndHessian(
1169 [[maybe_unused]] TMatrix const& F,
1170 [[maybe_unused]] typename TMatrix::ScalarType mu,
1171 [[maybe_unused]] typename TMatrix::ScalarType lambda,
1172 TMatrixGF& gF,
1173 TMatrixHF& HF) const
1174{
1175 static_assert(
1176 TMatrixGF::kRows == 9 and TMatrixGF::kCols == 1,
1177 "Grad w.r.t. F must have dimensions 9x1");
1178 static_assert(
1179 TMatrixHF::kRows == 9 and TMatrixHF::kCols == 9,
1180 "Hessian w.r.t. F must have dimensions 9x9");
1181 using ScalarType = typename TMatrix::ScalarType;
1182 ScalarType const a0 = F[4] * F[8];
1183 ScalarType const a1 = F[5] * F[7];
1184 ScalarType const a2 = 2 * a0 - 2 * a1;
1185 ScalarType const a3 = F[3] * F[8];
1186 ScalarType const a4 = F[4] * F[6];
1187 ScalarType const a5 = lambda * (-a1 * F[0] - a3 * F[1] - a4 * F[2] + F[0] * F[4] * F[8] +
1188 F[1] * F[5] * F[6] + F[2] * F[3] * F[7] - 1 - mu / lambda);
1189 ScalarType const a6 = (1.0 / 2.0) * a5;
1190 ScalarType const a7 = -2 * a3 + 2 * F[5] * F[6];
1191 ScalarType const a8 = F[3] * F[7];
1192 ScalarType const a9 = -2 * a4 + 2 * a8;
1193 ScalarType const a10 = F[1] * F[8];
1194 ScalarType const a11 = -2 * a10 + 2 * F[2] * F[7];
1195 ScalarType const a12 = F[0] * F[8];
1196 ScalarType const a13 = F[2] * F[6];
1197 ScalarType const a14 = 2 * a12 - 2 * a13;
1198 ScalarType const a15 = F[0] * F[7];
1199 ScalarType const a16 = -2 * a15 + 2 * F[1] * F[6];
1200 ScalarType const a17 = F[1] * F[5];
1201 ScalarType const a18 = F[2] * F[4];
1202 ScalarType const a19 = 2 * a17 - 2 * a18;
1203 ScalarType const a20 = F[0] * F[5];
1204 ScalarType const a21 = -2 * a20 + 2 * F[2] * F[3];
1205 ScalarType const a22 = F[0] * F[4];
1206 ScalarType const a23 = F[1] * F[3];
1207 ScalarType const a24 = 2 * a22 - 2 * a23;
1208 ScalarType const a25 = (1.0 / 2.0) * lambda;
1209 ScalarType const a26 = a25 * (a0 - a1);
1210 ScalarType const a27 = a5 * F[8];
1211 ScalarType const a28 = a5 * F[7];
1212 ScalarType const a29 = -a28;
1213 ScalarType const a30 = a5 * F[5];
1214 ScalarType const a31 = -a30;
1215 ScalarType const a32 = a5 * F[4];
1216 ScalarType const a33 = a25 * (-a3 + F[5] * F[6]);
1217 ScalarType const a34 = -a27;
1218 ScalarType const a35 = a5 * F[6];
1219 ScalarType const a36 = a5 * F[3];
1220 ScalarType const a37 = -a36;
1221 ScalarType const a38 = a25 * (-a4 + a8);
1222 ScalarType const a39 = -a35;
1223 ScalarType const a40 = -a32;
1224 ScalarType const a41 = a25 * (-a10 + F[2] * F[7]);
1225 ScalarType const a42 = a5 * F[2];
1226 ScalarType const a43 = a5 * F[1];
1227 ScalarType const a44 = -a43;
1228 ScalarType const a45 = a25 * (a12 - a13);
1229 ScalarType const a46 = -a42;
1230 ScalarType const a47 = a5 * F[0];
1231 ScalarType const a48 = a25 * (-a15 + F[1] * F[6]);
1232 ScalarType const a49 = -a47;
1233 ScalarType const a50 = a25 * (a17 - a18);
1234 ScalarType const a51 = a25 * (-a20 + F[2] * F[3]);
1235 ScalarType const a52 = a25 * (a22 - a23);
1236 gF[0] = a2 * a6 + mu * F[0];
1237 gF[1] = a6 * a7 + mu * F[1];
1238 gF[2] = a6 * a9 + mu * F[2];
1239 gF[3] = a11 * a6 + mu * F[3];
1240 gF[4] = a14 * a6 + mu * F[4];
1241 gF[5] = a16 * a6 + mu * F[5];
1242 gF[6] = a19 * a6 + mu * F[6];
1243 gF[7] = a21 * a6 + mu * F[7];
1244 gF[8] = a24 * a6 + mu * F[8];
1245 HF[0] = a2 * a26 + mu;
1246 HF[1] = a26 * a7;
1247 HF[2] = a26 * a9;
1248 HF[3] = a11 * a26;
1249 HF[4] = a14 * a26 + a27;
1250 HF[5] = a16 * a26 + a29;
1251 HF[6] = a19 * a26;
1252 HF[7] = a21 * a26 + a31;
1253 HF[8] = a24 * a26 + a32;
1254 HF[9] = a2 * a33;
1255 HF[10] = a33 * a7 + mu;
1256 HF[11] = a33 * a9;
1257 HF[12] = a11 * a33 + a34;
1258 HF[13] = a14 * a33;
1259 HF[14] = a16 * a33 + a35;
1260 HF[15] = a19 * a33 + a30;
1261 HF[16] = a21 * a33;
1262 HF[17] = a24 * a33 + a37;
1263 HF[18] = a2 * a38;
1264 HF[19] = a38 * a7;
1265 HF[20] = a38 * a9 + mu;
1266 HF[21] = a11 * a38 + a28;
1267 HF[22] = a14 * a38 + a39;
1268 HF[23] = a16 * a38;
1269 HF[24] = a19 * a38 + a40;
1270 HF[25] = a21 * a38 + a36;
1271 HF[26] = a24 * a38;
1272 HF[27] = a2 * a41;
1273 HF[28] = a34 + a41 * a7;
1274 HF[29] = a28 + a41 * a9;
1275 HF[30] = a11 * a41 + mu;
1276 HF[31] = a14 * a41;
1277 HF[32] = a16 * a41;
1278 HF[33] = a19 * a41;
1279 HF[34] = a21 * a41 + a42;
1280 HF[35] = a24 * a41 + a44;
1281 HF[36] = a2 * a45 + a27;
1282 HF[37] = a45 * a7;
1283 HF[38] = a39 + a45 * a9;
1284 HF[39] = a11 * a45;
1285 HF[40] = a14 * a45 + mu;
1286 HF[41] = a16 * a45;
1287 HF[42] = a19 * a45 + a46;
1288 HF[43] = a21 * a45;
1289 HF[44] = a24 * a45 + a47;
1290 HF[45] = a2 * a48 + a29;
1291 HF[46] = a35 + a48 * a7;
1292 HF[47] = a48 * a9;
1293 HF[48] = a11 * a48;
1294 HF[49] = a14 * a48;
1295 HF[50] = a16 * a48 + mu;
1296 HF[51] = a19 * a48 + a43;
1297 HF[52] = a21 * a48 + a49;
1298 HF[53] = a24 * a48;
1299 HF[54] = a2 * a50;
1300 HF[55] = a30 + a50 * a7;
1301 HF[56] = a40 + a50 * a9;
1302 HF[57] = a11 * a50;
1303 HF[58] = a14 * a50 + a46;
1304 HF[59] = a16 * a50 + a43;
1305 HF[60] = a19 * a50 + mu;
1306 HF[61] = a21 * a50;
1307 HF[62] = a24 * a50;
1308 HF[63] = a2 * a51 + a31;
1309 HF[64] = a51 * a7;
1310 HF[65] = a36 + a51 * a9;
1311 HF[66] = a11 * a51 + a42;
1312 HF[67] = a14 * a51;
1313 HF[68] = a16 * a51 + a49;
1314 HF[69] = a19 * a51;
1315 HF[70] = a21 * a51 + mu;
1316 HF[71] = a24 * a51;
1317 HF[72] = a2 * a52 + a32;
1318 HF[73] = a37 + a52 * a7;
1319 HF[74] = a52 * a9;
1320 HF[75] = a11 * a52 + a44;
1321 HF[76] = a14 * a52 + a47;
1322 HF[77] = a16 * a52;
1323 HF[78] = a19 * a52;
1324 HF[79] = a21 * a52;
1325 HF[80] = a24 * a52 + mu;
1326}
1327
1328} // namespace physics
1329} // namespace pbat
1330
1331#endif // PBAT_PHYSICS_STABLENEOHOOKEANENERGY_H
pbat::math::linalg::mini::SMatrix
Definition Matrix.h:121
pbat
The main namespace of the library.
Definition Aliases.h:15
pbat::physics::StableNeoHookeanEnergy< 1 >::evalWithGradAndHessian
PBAT_HOST_DEVICE TMatrix::ScalarType evalWithGradAndHessian(TMatrix const &F, typename TMatrix::ScalarType mu, typename TMatrix::ScalarType lambda, TMatrixGF &gF, TMatrixHF &HF) const
Evaluate the elastic energy with its gradient and hessian.
Definition StableNeoHookeanEnergy.h:229
pbat::physics::StableNeoHookeanEnergy< 1 >::hessian
PBAT_HOST_DEVICE SMatrix< typename TMatrix::ScalarType, 1, 1 > hessian(TMatrix const &F, typename TMatrix::ScalarType mu, typename TMatrix::ScalarType lambda) const
Evaluate the elastic energy hessian.
pbat::physics::StableNeoHookeanEnergy< 1 >::SMatrix
pbat::math::linalg::mini::SMatrix< TScalar, M, N > SMatrix
Scalar matrix type.
Definition StableNeoHookeanEnergy.h:35
pbat::physics::StableNeoHookeanEnergy< 1 >::kDims
static auto constexpr kDims
Dimension of the space.
Definition StableNeoHookeanEnergy.h:40
pbat::physics::StableNeoHookeanEnergy< 1 >::SVector
pbat::math::linalg::mini::SVector< TScalar, M > SVector
Scalar vector type.
Definition StableNeoHookeanEnergy.h:38
pbat::physics::StableNeoHookeanEnergy< 1 >::grad
PBAT_HOST_DEVICE SVector< typename TMatrix::ScalarType, 1 > grad(TMatrix const &F, typename TMatrix::ScalarType mu, typename TMatrix::ScalarType lambda) const
Evaluate the elastic energy gradient.
pbat::physics::StableNeoHookeanEnergy< 1 >::gradAndHessian
PBAT_HOST_DEVICE void gradAndHessian(TMatrix const &F, typename TMatrix::ScalarType mu, typename TMatrix::ScalarType lambda, TMatrixGF &gF, TMatrixHF &HF) const
Evaluate the elastic energy gradient and hessian.
Definition StableNeoHookeanEnergy.h:256
pbat::physics::StableNeoHookeanEnergy< 1 >::eval
PBAT_HOST_DEVICE TMatrix::ScalarType eval(TMatrix const &F, typename TMatrix::ScalarType mu, typename TMatrix::ScalarType lambda) const
Evaluate the elastic energy.
Definition StableNeoHookeanEnergy.h:148
pbat::physics::StableNeoHookeanEnergy< 1 >::evalWithGrad
PBAT_HOST_DEVICE TMatrix::ScalarType evalWithGrad(TMatrix const &F, typename TMatrix::ScalarType mu, typename TMatrix::ScalarType lambda, TMatrixGF &gF) const
Evaluate the elastic energy and its gradient.
Definition StableNeoHookeanEnergy.h:207
pbat::physics::StableNeoHookeanEnergy< 2 >::evalWithGradAndHessian
PBAT_HOST_DEVICE TMatrix::ScalarType evalWithGradAndHessian(TMatrix const &F, typename TMatrix::ScalarType mu, typename TMatrix::ScalarType lambda, TMatrixGF &gF, TMatrixHF &HF) const
Evaluate the elastic energy with its gradient and hessian.
Definition StableNeoHookeanEnergy.h:513
pbat::physics::StableNeoHookeanEnergy< 2 >::SMatrix
pbat::math::linalg::mini::SMatrix< TScalar, M, N > SMatrix
Scalar matrix type.
Definition StableNeoHookeanEnergy.h:283
pbat::physics::StableNeoHookeanEnergy< 2 >::gradAndHessian
PBAT_HOST_DEVICE void gradAndHessian(TMatrix const &F, typename TMatrix::ScalarType mu, typename TMatrix::ScalarType lambda, TMatrixGF &gF, TMatrixHF &HF) const
Evaluate the elastic energy gradient and hessian.
Definition StableNeoHookeanEnergy.h:570
pbat::physics::StableNeoHookeanEnergy< 2 >::grad
PBAT_HOST_DEVICE SVector< typename TMatrix::ScalarType, 4 > grad(TMatrix const &F, typename TMatrix::ScalarType mu, typename TMatrix::ScalarType lambda) const
Evaluate the elastic energy gradient.
pbat::physics::StableNeoHookeanEnergy< 2 >::hessian
PBAT_HOST_DEVICE SMatrix< typename TMatrix::ScalarType, 4, 4 > hessian(TMatrix const &F, typename TMatrix::ScalarType mu, typename TMatrix::ScalarType lambda) const
Evaluate the elastic energy hessian.
pbat::physics::StableNeoHookeanEnergy< 2 >::eval
PBAT_HOST_DEVICE TMatrix::ScalarType eval(TMatrix const &F, typename TMatrix::ScalarType mu, typename TMatrix::ScalarType lambda) const
Evaluate the elastic energy.
Definition StableNeoHookeanEnergy.h:396
pbat::physics::StableNeoHookeanEnergy< 2 >::SVector
pbat::math::linalg::mini::SVector< TScalar, M > SVector
Scalar vector type.
Definition StableNeoHookeanEnergy.h:286
pbat::physics::StableNeoHookeanEnergy< 2 >::kDims
static auto constexpr kDims
Dimension of the space.
Definition StableNeoHookeanEnergy.h:288
pbat::physics::StableNeoHookeanEnergy< 2 >::evalWithGrad
PBAT_HOST_DEVICE TMatrix::ScalarType evalWithGrad(TMatrix const &F, typename TMatrix::ScalarType mu, typename TMatrix::ScalarType lambda, TMatrixGF &gF) const
Evaluate the elastic energy and its gradient.
Definition StableNeoHookeanEnergy.h:486
pbat::physics::StableNeoHookeanEnergy< 3 >::kDims
static auto constexpr kDims
Dimension of the space.
Definition StableNeoHookeanEnergy.h:630
pbat::physics::StableNeoHookeanEnergy< 3 >::evalWithGrad
PBAT_HOST_DEVICE TMatrix::ScalarType evalWithGrad(TMatrix const &F, typename TMatrix::ScalarType mu, typename TMatrix::ScalarType lambda, TMatrixGF &gF) const
Evaluate the elastic energy and its gradient.
Definition StableNeoHookeanEnergy.h:955
pbat::physics::StableNeoHookeanEnergy< 3 >::evalWithGradAndHessian
PBAT_HOST_DEVICE TMatrix::ScalarType evalWithGradAndHessian(TMatrix const &F, typename TMatrix::ScalarType mu, typename TMatrix::ScalarType lambda, TMatrixGF &gF, TMatrixHF &HF) const
Evaluate the elastic energy with its gradient and hessian.
Definition StableNeoHookeanEnergy.h:997
pbat::physics::StableNeoHookeanEnergy< 3 >::SMatrix
pbat::math::linalg::mini::SMatrix< TScalar, M, N > SMatrix
Scalar matrix type.
Definition StableNeoHookeanEnergy.h:625
pbat::physics::StableNeoHookeanEnergy< 3 >::hessian
PBAT_HOST_DEVICE SMatrix< typename TMatrix::ScalarType, 9, 9 > hessian(TMatrix const &F, typename TMatrix::ScalarType mu, typename TMatrix::ScalarType lambda) const
Evaluate the elastic energy hessian.
pbat::physics::StableNeoHookeanEnergy< 3 >::eval
PBAT_HOST_DEVICE TMatrix::ScalarType eval(TMatrix const &F, typename TMatrix::ScalarType mu, typename TMatrix::ScalarType lambda) const
Evaluate the elastic energy.
Definition StableNeoHookeanEnergy.h:738
pbat::physics::StableNeoHookeanEnergy< 3 >::SVector
pbat::math::linalg::mini::SVector< TScalar, M > SVector
Scalar vector type.
Definition StableNeoHookeanEnergy.h:628
pbat::physics::StableNeoHookeanEnergy< 3 >::grad
PBAT_HOST_DEVICE SVector< typename TMatrix::ScalarType, 9 > grad(TMatrix const &F, typename TMatrix::ScalarType mu, typename TMatrix::ScalarType lambda) const
Evaluate the elastic energy gradient.
pbat::physics::StableNeoHookeanEnergy< 3 >::gradAndHessian
PBAT_HOST_DEVICE void gradAndHessian(TMatrix const &F, typename TMatrix::ScalarType mu, typename TMatrix::ScalarType lambda, TMatrixGF &gF, TMatrixHF &HF) const
Evaluate the elastic energy gradient and hessian.
Definition StableNeoHookeanEnergy.h:1168
pbat::physics::StableNeoHookeanEnergy
Definition StableNeoHookeanEnergy.h:23