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
SaintVenantKirchhoffEnergy.h
Go to the documentation of this file.
1
9
10#ifndef PBAT_PHYSICS_SAINTVENANTKIRCHHOFFENERGY_H
11#define PBAT_PHYSICS_SAINTVENANTKIRCHHOFFENERGY_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 SaintVenantKirchhoffEnergy;
24
30template <>
31struct SaintVenantKirchhoffEnergy<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 SaintVenantKirchhoffEnergy<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 ScalarType const a0 =
156 (((1.0 / 2.0) * ((F[0]) * (F[0])) - 1.0 / 2.0) *
157 ((1.0 / 2.0) * ((F[0]) * (F[0])) - 1.0 / 2.0));
158 psi = (1.0 / 2.0) * a0 * lambda + a0 * mu;
159 return psi;
160}
161
171template <math::linalg::mini::CReadableVectorizedMatrix TMatrix>
172PBAT_HOST_DEVICE SaintVenantKirchhoffEnergy<1>::SVector<typename TMatrix::ScalarType, 1>
173SaintVenantKirchhoffEnergy<1>::grad(
174 [[maybe_unused]] TMatrix const& F,
175 [[maybe_unused]] typename TMatrix::ScalarType mu,
176 [[maybe_unused]] typename TMatrix::ScalarType lambda) const
177{
178 using ScalarType = typename TMatrix::ScalarType;
179 SVector<ScalarType, 1> G;
180 ScalarType const a0 = ((1.0 / 2.0) * ((F[0]) * (F[0])) - 1.0 / 2.0) * F[0];
181 G[0] = a0 * lambda + 2 * a0 * mu;
182 return G;
183}
184
194template <math::linalg::mini::CReadableVectorizedMatrix TMatrix>
195PBAT_HOST_DEVICE SaintVenantKirchhoffEnergy<1>::SMatrix<typename TMatrix::ScalarType, 1, 1>
196SaintVenantKirchhoffEnergy<1>::hessian(
197 [[maybe_unused]] TMatrix const& F,
198 [[maybe_unused]] typename TMatrix::ScalarType mu,
199 [[maybe_unused]] typename TMatrix::ScalarType lambda) const
200{
201 using ScalarType = typename TMatrix::ScalarType;
202 SMatrix<ScalarType, 1, 1> H;
203 ScalarType const a0 = ((F[0]) * (F[0]));
204 ScalarType const a1 = 2 * mu;
205 ScalarType const a2 = (1.0 / 2.0) * a0 - 1.0 / 2.0;
206 H[0] = a0 * a1 + a0 * lambda + a1 * a2 + a2 * lambda;
207 return H;
208}
209
210template <
211 math::linalg::mini::CReadableVectorizedMatrix TMatrix,
212 math::linalg::mini::CWriteableVectorizedMatrix TMatrixGF>
213PBAT_HOST_DEVICE typename TMatrix::ScalarType SaintVenantKirchhoffEnergy<1>::evalWithGrad(
214 [[maybe_unused]] TMatrix const& F,
215 [[maybe_unused]] typename TMatrix::ScalarType mu,
216 [[maybe_unused]] typename TMatrix::ScalarType lambda,
217 TMatrixGF& gF) const
218{
219 static_assert(
220 TMatrixGF::kRows == 1 and TMatrixGF::kCols == 1,
221 "Grad w.r.t. F must have dimensions 1x1");
222 using ScalarType = typename TMatrix::ScalarType;
223 ScalarType psi;
224 ScalarType const a0 = (1.0 / 2.0) * ((F[0]) * (F[0])) - 1.0 / 2.0;
225 ScalarType const a1 = ((a0) * (a0));
226 ScalarType const a2 = a0 * F[0];
227 psi = (1.0 / 2.0) * a1 * lambda + a1 * mu;
228 gF[0] = a2 * lambda + 2 * a2 * mu;
229 return psi;
230}
231
232template <
233 math::linalg::mini::CReadableVectorizedMatrix TMatrix,
234 math::linalg::mini::CWriteableVectorizedMatrix TMatrixGF,
235 math::linalg::mini::CWriteableVectorizedMatrix TMatrixHF>
236PBAT_HOST_DEVICE typename TMatrix::ScalarType SaintVenantKirchhoffEnergy<1>::evalWithGradAndHessian(
237 [[maybe_unused]] TMatrix const& F,
238 [[maybe_unused]] typename TMatrix::ScalarType mu,
239 [[maybe_unused]] typename TMatrix::ScalarType lambda,
240 TMatrixGF& gF,
241 TMatrixHF& HF) const
242{
243 static_assert(
244 TMatrixGF::kRows == 1 and TMatrixGF::kCols == 1,
245 "Grad w.r.t. F must have dimensions 1x1");
246 static_assert(
247 TMatrixHF::kRows == 1 and TMatrixHF::kCols == 1,
248 "Hessian w.r.t. F must have dimensions 1x1");
249 using ScalarType = typename TMatrix::ScalarType;
250 ScalarType psi;
251 ScalarType const a0 = ((F[0]) * (F[0]));
252 ScalarType const a1 = (1.0 / 2.0) * a0 - 1.0 / 2.0;
253 ScalarType const a2 = ((a1) * (a1));
254 ScalarType const a3 = a1 * lambda;
255 ScalarType const a4 = 2 * mu;
256 ScalarType const a5 = a1 * a4;
257 psi = (1.0 / 2.0) * a2 * lambda + a2 * mu;
258 gF[0] = a3 * F[0] + a5 * F[0];
259 HF[0] = a0 * a4 + a0 * lambda + a3 + a5;
260 return psi;
261}
262
263template <
264 math::linalg::mini::CReadableVectorizedMatrix TMatrix,
265 math::linalg::mini::CWriteableVectorizedMatrix TMatrixGF,
266 math::linalg::mini::CWriteableVectorizedMatrix TMatrixHF>
267PBAT_HOST_DEVICE void SaintVenantKirchhoffEnergy<1>::gradAndHessian(
268 [[maybe_unused]] TMatrix const& F,
269 [[maybe_unused]] typename TMatrix::ScalarType mu,
270 [[maybe_unused]] typename TMatrix::ScalarType lambda,
271 TMatrixGF& gF,
272 TMatrixHF& HF) const
273{
274 static_assert(
275 TMatrixGF::kRows == 1 and TMatrixGF::kCols == 1,
276 "Grad w.r.t. F must have dimensions 1x1");
277 static_assert(
278 TMatrixHF::kRows == 1 and TMatrixHF::kCols == 1,
279 "Hessian w.r.t. F must have dimensions 1x1");
280 using ScalarType = typename TMatrix::ScalarType;
281 ScalarType const a0 = ((F[0]) * (F[0]));
282 ScalarType const a1 = (1.0 / 2.0) * a0 - 1.0 / 2.0;
283 ScalarType const a2 = a1 * lambda;
284 ScalarType const a3 = 2 * mu;
285 ScalarType const a4 = a1 * a3;
286 gF[0] = a2 * F[0] + a4 * F[0];
287 HF[0] = a0 * a3 + a0 * lambda + a2 + a4;
288}
289
295template <>
296struct SaintVenantKirchhoffEnergy<2>
297{
298 public:
299 template <class TScalar, int M, int N>
300 using SMatrix = pbat::math::linalg::mini::SMatrix<TScalar, M, N>;
301
302 template <class TScalar, int M>
303 using SVector = pbat::math::linalg::mini::SVector<TScalar, M>;
304
305 static auto constexpr kDims = 2;
306
316 template <math::linalg::mini::CReadableVectorizedMatrix TMatrix>
317 PBAT_HOST_DEVICE typename TMatrix::ScalarType
318 eval(TMatrix const& F, typename TMatrix::ScalarType mu, typename TMatrix::ScalarType lambda)
319 const;
320
330 template <math::linalg::mini::CReadableVectorizedMatrix TMatrix>
331 PBAT_HOST_DEVICE SVector<typename TMatrix::ScalarType, 4>
332 grad(TMatrix const& F, typename TMatrix::ScalarType mu, typename TMatrix::ScalarType lambda)
333 const;
334
344 template <math::linalg::mini::CReadableVectorizedMatrix TMatrix>
345 PBAT_HOST_DEVICE SMatrix<typename TMatrix::ScalarType, 4, 4>
346 hessian(TMatrix const& F, typename TMatrix::ScalarType mu, typename TMatrix::ScalarType lambda)
347 const;
348
359 template <
360 math::linalg::mini::CReadableVectorizedMatrix TMatrix,
361 math::linalg::mini::CWriteableVectorizedMatrix TMatrixGF>
362 PBAT_HOST_DEVICE typename TMatrix::ScalarType evalWithGrad(
363 TMatrix const& F,
364 typename TMatrix::ScalarType mu,
365 typename TMatrix::ScalarType lambda,
366 TMatrixGF& gF) const;
367
379 template <
380 math::linalg::mini::CReadableVectorizedMatrix TMatrix,
381 math::linalg::mini::CWriteableVectorizedMatrix TMatrixGF,
382 math::linalg::mini::CWriteableVectorizedMatrix TMatrixHF>
383 PBAT_HOST_DEVICE typename TMatrix::ScalarType evalWithGradAndHessian(
384 TMatrix const& F,
385 typename TMatrix::ScalarType mu,
386 typename TMatrix::ScalarType lambda,
387 TMatrixGF& gF,
388 TMatrixHF& HF) const;
389
400 template <
401 math::linalg::mini::CReadableVectorizedMatrix TMatrix,
402 math::linalg::mini::CWriteableVectorizedMatrix TMatrixGF,
403 math::linalg::mini::CWriteableVectorizedMatrix TMatrixHF>
404 PBAT_HOST_DEVICE void gradAndHessian(
405 TMatrix const& F,
406 typename TMatrix::ScalarType mu,
407 typename TMatrix::ScalarType lambda,
408 TMatrixGF& gF,
409 TMatrixHF& HF) const;
410};
411
412template <math::linalg::mini::CReadableVectorizedMatrix TMatrix>
413PBAT_HOST_DEVICE typename TMatrix::ScalarType SaintVenantKirchhoffEnergy<2>::eval(
414 [[maybe_unused]] TMatrix const& F,
415 [[maybe_unused]] typename TMatrix::ScalarType mu,
416 [[maybe_unused]] typename TMatrix::ScalarType lambda) const
417{
418 using ScalarType = typename TMatrix::ScalarType;
419 ScalarType psi;
420 ScalarType const a0 = (1.0 / 2.0) * ((F[0]) * (F[0])) + (1.0 / 2.0) * ((F[1]) * (F[1]));
421 ScalarType const a1 = (1.0 / 2.0) * ((F[2]) * (F[2])) + (1.0 / 2.0) * ((F[3]) * (F[3]));
422 psi = (1.0 / 2.0) * lambda * ((a0 + a1 - 1) * (a0 + a1 - 1)) +
423 mu * (((a0 - 1.0 / 2.0) * (a0 - 1.0 / 2.0)) + ((a1 - 1.0 / 2.0) * (a1 - 1.0 / 2.0)) +
424 2 * (((1.0 / 2.0) * F[0] * F[2] + (1.0 / 2.0) * F[1] * F[3]) *
425 ((1.0 / 2.0) * F[0] * F[2] + (1.0 / 2.0) * F[1] * F[3])));
426 return psi;
427}
428
438template <math::linalg::mini::CReadableVectorizedMatrix TMatrix>
439PBAT_HOST_DEVICE SaintVenantKirchhoffEnergy<2>::SVector<typename TMatrix::ScalarType, 4>
440SaintVenantKirchhoffEnergy<2>::grad(
441 [[maybe_unused]] TMatrix const& F,
442 [[maybe_unused]] typename TMatrix::ScalarType mu,
443 [[maybe_unused]] typename TMatrix::ScalarType lambda) const
444{
445 using ScalarType = typename TMatrix::ScalarType;
446 SVector<ScalarType, 4> G;
447 ScalarType const a0 = (1.0 / 2.0) * ((F[0]) * (F[0])) + (1.0 / 2.0) * ((F[1]) * (F[1]));
448 ScalarType const a1 = (1.0 / 2.0) * ((F[2]) * (F[2])) + (1.0 / 2.0) * ((F[3]) * (F[3]));
449 ScalarType const a2 = lambda * (a0 + a1 - 1);
450 ScalarType const a3 = 2 * a0 - 1;
451 ScalarType const a4 = F[0] * F[2] + F[1] * F[3];
452 ScalarType const a5 = 2 * a1 - 1;
453 G[0] = a2 * F[0] + mu * (a3 * F[0] + a4 * F[2]);
454 G[1] = a2 * F[1] + mu * (a3 * F[1] + a4 * F[3]);
455 G[2] = a2 * F[2] + mu * (a4 * F[0] + a5 * F[2]);
456 G[3] = a2 * F[3] + mu * (a4 * F[1] + a5 * F[3]);
457 return G;
458}
459
469template <math::linalg::mini::CReadableVectorizedMatrix TMatrix>
470PBAT_HOST_DEVICE SaintVenantKirchhoffEnergy<2>::SMatrix<typename TMatrix::ScalarType, 4, 4>
471SaintVenantKirchhoffEnergy<2>::hessian(
472 [[maybe_unused]] TMatrix const& F,
473 [[maybe_unused]] typename TMatrix::ScalarType mu,
474 [[maybe_unused]] typename TMatrix::ScalarType lambda) const
475{
476 using ScalarType = typename TMatrix::ScalarType;
477 SMatrix<ScalarType, 4, 4> H;
478 ScalarType const a0 = ((F[0]) * (F[0]));
479 ScalarType const a1 = ((F[1]) * (F[1]));
480 ScalarType const a2 = ((F[2]) * (F[2]));
481 ScalarType const a3 = a1 + a2 - 1;
482 ScalarType const a4 = ((F[3]) * (F[3]));
483 ScalarType const a5 =
484 lambda * ((1.0 / 2.0) * a0 + (1.0 / 2.0) * a1 + (1.0 / 2.0) * a2 + (1.0 / 2.0) * a4 - 1);
485 ScalarType const a6 = F[0] * F[1];
486 ScalarType const a7 = F[2] * F[3];
487 ScalarType const a8 = a6 * lambda + mu * (2 * a6 + a7);
488 ScalarType const a9 = F[0] * F[2];
489 ScalarType const a10 = F[1] * F[3];
490 ScalarType const a11 = a9 * lambda + mu * (a10 + 2 * a9);
491 ScalarType const a12 = F[0] * F[3];
492 ScalarType const a13 = F[1] * F[2];
493 ScalarType const a14 = a12 * lambda + a13 * mu;
494 ScalarType const a15 = a0 + a4 - 1;
495 ScalarType const a16 = a12 * mu + a13 * lambda;
496 ScalarType const a17 = a10 * lambda + mu * (2 * a10 + a9);
497 ScalarType const a18 = a7 * lambda + mu * (a6 + 2 * a7);
498 H[0] = a0 * lambda + a5 + mu * (3 * a0 + a3);
499 H[1] = a8;
500 H[2] = a11;
501 H[3] = a14;
502 H[4] = a8;
503 H[5] = a1 * lambda + a5 + mu * (3 * a1 + a15);
504 H[6] = a16;
505 H[7] = a17;
506 H[8] = a11;
507 H[9] = a16;
508 H[10] = a2 * lambda + a5 + mu * (a15 + 3 * a2);
509 H[11] = a18;
510 H[12] = a14;
511 H[13] = a17;
512 H[14] = a18;
513 H[15] = a4 * lambda + a5 + mu * (a3 + 3 * a4);
514 return H;
515}
516
517template <
518 math::linalg::mini::CReadableVectorizedMatrix TMatrix,
519 math::linalg::mini::CWriteableVectorizedMatrix TMatrixGF>
520PBAT_HOST_DEVICE typename TMatrix::ScalarType SaintVenantKirchhoffEnergy<2>::evalWithGrad(
521 [[maybe_unused]] TMatrix const& F,
522 [[maybe_unused]] typename TMatrix::ScalarType mu,
523 [[maybe_unused]] typename TMatrix::ScalarType lambda,
524 TMatrixGF& gF) const
525{
526 static_assert(
527 TMatrixGF::kRows == 4 and TMatrixGF::kCols == 1,
528 "Grad w.r.t. F must have dimensions 4x1");
529 using ScalarType = typename TMatrix::ScalarType;
530 ScalarType psi;
531 ScalarType const a0 = (1.0 / 2.0) * ((F[0]) * (F[0])) + (1.0 / 2.0) * ((F[1]) * (F[1]));
532 ScalarType const a1 = (1.0 / 2.0) * ((F[2]) * (F[2])) + (1.0 / 2.0) * ((F[3]) * (F[3]));
533 ScalarType const a2 = a0 + a1 - 1;
534 ScalarType const a3 = a0 - 1.0 / 2.0;
535 ScalarType const a4 = a1 - 1.0 / 2.0;
536 ScalarType const a5 = (1.0 / 2.0) * F[0] * F[2] + (1.0 / 2.0) * F[1] * F[3];
537 ScalarType const a6 = a2 * lambda;
538 ScalarType const a7 = 2 * a3;
539 ScalarType const a8 = 2 * a5;
540 ScalarType const a9 = 2 * a4;
541 psi = (1.0 / 2.0) * ((a2) * (a2)) * lambda +
542 mu * (((a3) * (a3)) + ((a4) * (a4)) + 2 * ((a5) * (a5)));
543 gF[0] = a6 * F[0] + mu * (a7 * F[0] + a8 * F[2]);
544 gF[1] = a6 * F[1] + mu * (a7 * F[1] + a8 * F[3]);
545 gF[2] = a6 * F[2] + mu * (a8 * F[0] + a9 * F[2]);
546 gF[3] = a6 * F[3] + mu * (a8 * F[1] + a9 * F[3]);
547 return psi;
548}
549
550template <
551 math::linalg::mini::CReadableVectorizedMatrix TMatrix,
552 math::linalg::mini::CWriteableVectorizedMatrix TMatrixGF,
553 math::linalg::mini::CWriteableVectorizedMatrix TMatrixHF>
554PBAT_HOST_DEVICE typename TMatrix::ScalarType SaintVenantKirchhoffEnergy<2>::evalWithGradAndHessian(
555 [[maybe_unused]] TMatrix const& F,
556 [[maybe_unused]] typename TMatrix::ScalarType mu,
557 [[maybe_unused]] typename TMatrix::ScalarType lambda,
558 TMatrixGF& gF,
559 TMatrixHF& HF) const
560{
561 static_assert(
562 TMatrixGF::kRows == 4 and TMatrixGF::kCols == 1,
563 "Grad w.r.t. F must have dimensions 4x1");
564 static_assert(
565 TMatrixHF::kRows == 4 and TMatrixHF::kCols == 4,
566 "Hessian w.r.t. F must have dimensions 4x4");
567 using ScalarType = typename TMatrix::ScalarType;
568 ScalarType psi;
569 ScalarType const a0 = ((F[0]) * (F[0]));
570 ScalarType const a1 = ((F[1]) * (F[1]));
571 ScalarType const a2 = (1.0 / 2.0) * a0 + (1.0 / 2.0) * a1;
572 ScalarType const a3 = ((F[2]) * (F[2]));
573 ScalarType const a4 = ((F[3]) * (F[3]));
574 ScalarType const a5 = (1.0 / 2.0) * a3 + (1.0 / 2.0) * a4;
575 ScalarType const a6 = a2 + a5 - 1;
576 ScalarType const a7 = a2 - 1.0 / 2.0;
577 ScalarType const a8 = a5 - 1.0 / 2.0;
578 ScalarType const a9 = F[0] * F[2];
579 ScalarType const a10 = F[1] * F[3];
580 ScalarType const a11 = (1.0 / 2.0) * a10 + (1.0 / 2.0) * a9;
581 ScalarType const a12 = a6 * lambda;
582 ScalarType const a13 = 2 * a7;
583 ScalarType const a14 = 2 * a11;
584 ScalarType const a15 = 2 * a8;
585 ScalarType const a16 = a1 + a3 - 1;
586 ScalarType const a17 = F[0] * F[1];
587 ScalarType const a18 = F[2] * F[3];
588 ScalarType const a19 = a17 * lambda + mu * (2 * a17 + a18);
589 ScalarType const a20 = a9 * lambda + mu * (a10 + 2 * a9);
590 ScalarType const a21 = F[0] * F[3];
591 ScalarType const a22 = F[1] * F[2];
592 ScalarType const a23 = a21 * lambda + a22 * mu;
593 ScalarType const a24 = a0 + a4 - 1;
594 ScalarType const a25 = a21 * mu + a22 * lambda;
595 ScalarType const a26 = a10 * lambda + mu * (2 * a10 + a9);
596 ScalarType const a27 = a18 * lambda + mu * (a17 + 2 * a18);
597 psi = (1.0 / 2.0) * ((a6) * (a6)) * lambda +
598 mu * (2 * ((a11) * (a11)) + ((a7) * (a7)) + ((a8) * (a8)));
599 gF[0] = a12 * F[0] + mu * (a13 * F[0] + a14 * F[2]);
600 gF[1] = a12 * F[1] + mu * (a13 * F[1] + a14 * F[3]);
601 gF[2] = a12 * F[2] + mu * (a14 * F[0] + a15 * F[2]);
602 gF[3] = a12 * F[3] + mu * (a14 * F[1] + a15 * F[3]);
603 HF[0] = a0 * lambda + a12 + mu * (3 * a0 + a16);
604 HF[1] = a19;
605 HF[2] = a20;
606 HF[3] = a23;
607 HF[4] = a19;
608 HF[5] = a1 * lambda + a12 + mu * (3 * a1 + a24);
609 HF[6] = a25;
610 HF[7] = a26;
611 HF[8] = a20;
612 HF[9] = a25;
613 HF[10] = a12 + a3 * lambda + mu * (a24 + 3 * a3);
614 HF[11] = a27;
615 HF[12] = a23;
616 HF[13] = a26;
617 HF[14] = a27;
618 HF[15] = a12 + a4 * lambda + mu * (a16 + 3 * a4);
619 return psi;
620}
621
622template <
623 math::linalg::mini::CReadableVectorizedMatrix TMatrix,
624 math::linalg::mini::CWriteableVectorizedMatrix TMatrixGF,
625 math::linalg::mini::CWriteableVectorizedMatrix TMatrixHF>
626PBAT_HOST_DEVICE void SaintVenantKirchhoffEnergy<2>::gradAndHessian(
627 [[maybe_unused]] TMatrix const& F,
628 [[maybe_unused]] typename TMatrix::ScalarType mu,
629 [[maybe_unused]] typename TMatrix::ScalarType lambda,
630 TMatrixGF& gF,
631 TMatrixHF& HF) const
632{
633 static_assert(
634 TMatrixGF::kRows == 4 and TMatrixGF::kCols == 1,
635 "Grad w.r.t. F must have dimensions 4x1");
636 static_assert(
637 TMatrixHF::kRows == 4 and TMatrixHF::kCols == 4,
638 "Hessian w.r.t. F must have dimensions 4x4");
639 using ScalarType = typename TMatrix::ScalarType;
640 ScalarType const a0 = ((F[0]) * (F[0]));
641 ScalarType const a1 = ((F[1]) * (F[1]));
642 ScalarType const a2 = (1.0 / 2.0) * a0 + (1.0 / 2.0) * a1;
643 ScalarType const a3 = ((F[2]) * (F[2]));
644 ScalarType const a4 = ((F[3]) * (F[3]));
645 ScalarType const a5 = (1.0 / 2.0) * a3 + (1.0 / 2.0) * a4;
646 ScalarType const a6 = lambda * (a2 + a5 - 1);
647 ScalarType const a7 = 2 * a2 - 1;
648 ScalarType const a8 = F[0] * F[2];
649 ScalarType const a9 = F[1] * F[3];
650 ScalarType const a10 = a8 + a9;
651 ScalarType const a11 = 2 * a5 - 1;
652 ScalarType const a12 = a1 + a3 - 1;
653 ScalarType const a13 = F[0] * F[1];
654 ScalarType const a14 = F[2] * F[3];
655 ScalarType const a15 = a13 * lambda + mu * (2 * a13 + a14);
656 ScalarType const a16 = a8 * lambda + mu * (2 * a8 + a9);
657 ScalarType const a17 = F[0] * F[3];
658 ScalarType const a18 = F[1] * F[2];
659 ScalarType const a19 = a17 * lambda + a18 * mu;
660 ScalarType const a20 = a0 + a4 - 1;
661 ScalarType const a21 = a17 * mu + a18 * lambda;
662 ScalarType const a22 = a9 * lambda + mu * (a8 + 2 * a9);
663 ScalarType const a23 = a14 * lambda + mu * (a13 + 2 * a14);
664 gF[0] = a6 * F[0] + mu * (a10 * F[2] + a7 * F[0]);
665 gF[1] = a6 * F[1] + mu * (a10 * F[3] + a7 * F[1]);
666 gF[2] = a6 * F[2] + mu * (a10 * F[0] + a11 * F[2]);
667 gF[3] = a6 * F[3] + mu * (a10 * F[1] + a11 * F[3]);
668 HF[0] = a0 * lambda + a6 + mu * (3 * a0 + a12);
669 HF[1] = a15;
670 HF[2] = a16;
671 HF[3] = a19;
672 HF[4] = a15;
673 HF[5] = a1 * lambda + a6 + mu * (3 * a1 + a20);
674 HF[6] = a21;
675 HF[7] = a22;
676 HF[8] = a16;
677 HF[9] = a21;
678 HF[10] = a3 * lambda + a6 + mu * (a20 + 3 * a3);
679 HF[11] = a23;
680 HF[12] = a19;
681 HF[13] = a22;
682 HF[14] = a23;
683 HF[15] = a4 * lambda + a6 + mu * (a12 + 3 * a4);
684}
685
691template <>
692struct SaintVenantKirchhoffEnergy<3>
693{
694 public:
695 template <class TScalar, int M, int N>
696 using SMatrix = pbat::math::linalg::mini::SMatrix<TScalar, M, N>;
697
698 template <class TScalar, int M>
699 using SVector = pbat::math::linalg::mini::SVector<TScalar, M>;
700
701 static auto constexpr kDims = 3;
702
712 template <math::linalg::mini::CReadableVectorizedMatrix TMatrix>
713 PBAT_HOST_DEVICE typename TMatrix::ScalarType
714 eval(TMatrix const& F, typename TMatrix::ScalarType mu, typename TMatrix::ScalarType lambda)
715 const;
716
726 template <math::linalg::mini::CReadableVectorizedMatrix TMatrix>
727 PBAT_HOST_DEVICE SVector<typename TMatrix::ScalarType, 9>
728 grad(TMatrix const& F, typename TMatrix::ScalarType mu, typename TMatrix::ScalarType lambda)
729 const;
730
740 template <math::linalg::mini::CReadableVectorizedMatrix TMatrix>
741 PBAT_HOST_DEVICE SMatrix<typename TMatrix::ScalarType, 9, 9>
742 hessian(TMatrix const& F, typename TMatrix::ScalarType mu, typename TMatrix::ScalarType lambda)
743 const;
744
755 template <
756 math::linalg::mini::CReadableVectorizedMatrix TMatrix,
757 math::linalg::mini::CWriteableVectorizedMatrix TMatrixGF>
758 PBAT_HOST_DEVICE typename TMatrix::ScalarType evalWithGrad(
759 TMatrix const& F,
760 typename TMatrix::ScalarType mu,
761 typename TMatrix::ScalarType lambda,
762 TMatrixGF& gF) const;
763
775 template <
776 math::linalg::mini::CReadableVectorizedMatrix TMatrix,
777 math::linalg::mini::CWriteableVectorizedMatrix TMatrixGF,
778 math::linalg::mini::CWriteableVectorizedMatrix TMatrixHF>
779 PBAT_HOST_DEVICE typename TMatrix::ScalarType evalWithGradAndHessian(
780 TMatrix const& F,
781 typename TMatrix::ScalarType mu,
782 typename TMatrix::ScalarType lambda,
783 TMatrixGF& gF,
784 TMatrixHF& HF) const;
785
796 template <
797 math::linalg::mini::CReadableVectorizedMatrix TMatrix,
798 math::linalg::mini::CWriteableVectorizedMatrix TMatrixGF,
799 math::linalg::mini::CWriteableVectorizedMatrix TMatrixHF>
800 PBAT_HOST_DEVICE void gradAndHessian(
801 TMatrix const& F,
802 typename TMatrix::ScalarType mu,
803 typename TMatrix::ScalarType lambda,
804 TMatrixGF& gF,
805 TMatrixHF& HF) const;
806};
807
808template <math::linalg::mini::CReadableVectorizedMatrix TMatrix>
809PBAT_HOST_DEVICE typename TMatrix::ScalarType SaintVenantKirchhoffEnergy<3>::eval(
810 [[maybe_unused]] TMatrix const& F,
811 [[maybe_unused]] typename TMatrix::ScalarType mu,
812 [[maybe_unused]] typename TMatrix::ScalarType lambda) const
813{
814 using ScalarType = typename TMatrix::ScalarType;
815 ScalarType psi;
816 ScalarType const a0 = (1.0 / 2.0) * ((F[0]) * (F[0])) + (1.0 / 2.0) * ((F[1]) * (F[1])) +
817 (1.0 / 2.0) * ((F[2]) * (F[2]));
818 ScalarType const a1 = (1.0 / 2.0) * ((F[3]) * (F[3])) + (1.0 / 2.0) * ((F[4]) * (F[4])) +
819 (1.0 / 2.0) * ((F[5]) * (F[5]));
820 ScalarType const a2 = (1.0 / 2.0) * ((F[6]) * (F[6])) + (1.0 / 2.0) * ((F[7]) * (F[7])) +
821 (1.0 / 2.0) * ((F[8]) * (F[8]));
822 ScalarType const a3 = (1.0 / 2.0) * F[0];
823 ScalarType const a4 = (1.0 / 2.0) * F[1];
824 ScalarType const a5 = (1.0 / 2.0) * F[2];
825 psi = (1.0 / 2.0) * lambda * ((a0 + a1 + a2 - 3.0 / 2.0) * (a0 + a1 + a2 - 3.0 / 2.0)) +
826 mu * (((a0 - 1.0 / 2.0) * (a0 - 1.0 / 2.0)) + ((a1 - 1.0 / 2.0) * (a1 - 1.0 / 2.0)) +
827 ((a2 - 1.0 / 2.0) * (a2 - 1.0 / 2.0)) +
828 2 * ((a3 * F[3] + a4 * F[4] + a5 * F[5]) * (a3 * F[3] + a4 * F[4] + a5 * F[5])) +
829 2 * ((a3 * F[6] + a4 * F[7] + a5 * F[8]) * (a3 * F[6] + a4 * F[7] + a5 * F[8])) +
830 2 * (((1.0 / 2.0) * F[3] * F[6] + (1.0 / 2.0) * F[4] * F[7] +
831 (1.0 / 2.0) * F[5] * F[8]) *
832 ((1.0 / 2.0) * F[3] * F[6] + (1.0 / 2.0) * F[4] * F[7] +
833 (1.0 / 2.0) * F[5] * F[8])));
834 return psi;
835}
836
846template <math::linalg::mini::CReadableVectorizedMatrix TMatrix>
847PBAT_HOST_DEVICE SaintVenantKirchhoffEnergy<3>::SVector<typename TMatrix::ScalarType, 9>
848SaintVenantKirchhoffEnergy<3>::grad(
849 [[maybe_unused]] TMatrix const& F,
850 [[maybe_unused]] typename TMatrix::ScalarType mu,
851 [[maybe_unused]] typename TMatrix::ScalarType lambda) const
852{
853 using ScalarType = typename TMatrix::ScalarType;
854 SVector<ScalarType, 9> G;
855 ScalarType const a0 = (1.0 / 2.0) * ((F[0]) * (F[0])) + (1.0 / 2.0) * ((F[1]) * (F[1])) +
856 (1.0 / 2.0) * ((F[2]) * (F[2]));
857 ScalarType const a1 = (1.0 / 2.0) * ((F[3]) * (F[3])) + (1.0 / 2.0) * ((F[4]) * (F[4])) +
858 (1.0 / 2.0) * ((F[5]) * (F[5]));
859 ScalarType const a2 = (1.0 / 2.0) * ((F[6]) * (F[6])) + (1.0 / 2.0) * ((F[7]) * (F[7])) +
860 (1.0 / 2.0) * ((F[8]) * (F[8]));
861 ScalarType const a3 = lambda * (a0 + a1 + a2 - 3.0 / 2.0);
862 ScalarType const a4 = 2 * a0 - 1;
863 ScalarType const a5 = (1.0 / 2.0) * F[0];
864 ScalarType const a6 = (1.0 / 2.0) * F[1];
865 ScalarType const a7 = (1.0 / 2.0) * F[2];
866 ScalarType const a8 = 2 * a5 * F[3] + 2 * a6 * F[4] + 2 * a7 * F[5];
867 ScalarType const a9 = 2 * a5 * F[6] + 2 * a6 * F[7] + 2 * a7 * F[8];
868 ScalarType const a10 = 2 * a1 - 1;
869 ScalarType const a11 = F[3] * F[6] + F[4] * F[7] + F[5] * F[8];
870 ScalarType const a12 = 2 * a2 - 1;
871 G[0] = a3 * F[0] + mu * (a4 * F[0] + a8 * F[3] + a9 * F[6]);
872 G[1] = a3 * F[1] + mu * (a4 * F[1] + a8 * F[4] + a9 * F[7]);
873 G[2] = a3 * F[2] + mu * (a4 * F[2] + a8 * F[5] + a9 * F[8]);
874 G[3] = a3 * F[3] + mu * (a10 * F[3] + a11 * F[6] + a8 * F[0]);
875 G[4] = a3 * F[4] + mu * (a10 * F[4] + a11 * F[7] + a8 * F[1]);
876 G[5] = a3 * F[5] + mu * (a10 * F[5] + a11 * F[8] + a8 * F[2]);
877 G[6] = a3 * F[6] + mu * (a11 * F[3] + a12 * F[6] + a9 * F[0]);
878 G[7] = a3 * F[7] + mu * (a11 * F[4] + a12 * F[7] + a9 * F[1]);
879 G[8] = a3 * F[8] + mu * (a11 * F[5] + a12 * F[8] + a9 * F[2]);
880 return G;
881}
882
892template <math::linalg::mini::CReadableVectorizedMatrix TMatrix>
893PBAT_HOST_DEVICE SaintVenantKirchhoffEnergy<3>::SMatrix<typename TMatrix::ScalarType, 9, 9>
894SaintVenantKirchhoffEnergy<3>::hessian(
895 [[maybe_unused]] TMatrix const& F,
896 [[maybe_unused]] typename TMatrix::ScalarType mu,
897 [[maybe_unused]] typename TMatrix::ScalarType lambda) const
898{
899 using ScalarType = typename TMatrix::ScalarType;
900 SMatrix<ScalarType, 9, 9> H;
901 ScalarType const a0 = ((F[0]) * (F[0]));
902 ScalarType const a1 = ((F[1]) * (F[1]));
903 ScalarType const a2 = ((F[3]) * (F[3]));
904 ScalarType const a3 = a1 + a2;
905 ScalarType const a4 = ((F[6]) * (F[6]));
906 ScalarType const a5 = ((F[2]) * (F[2]));
907 ScalarType const a6 = a5 - 1;
908 ScalarType const a7 = a4 + a6;
909 ScalarType const a8 = ((F[4]) * (F[4]));
910 ScalarType const a9 = ((F[5]) * (F[5]));
911 ScalarType const a10 = ((F[7]) * (F[7]));
912 ScalarType const a11 = ((F[8]) * (F[8]));
913 ScalarType const a12 =
914 lambda * ((1.0 / 2.0) * a0 + (1.0 / 2.0) * a1 + (1.0 / 2.0) * a10 + (1.0 / 2.0) * a11 +
915 (1.0 / 2.0) * a2 + (1.0 / 2.0) * a4 + (1.0 / 2.0) * a5 + (1.0 / 2.0) * a8 +
916 (1.0 / 2.0) * a9 - 3.0 / 2.0);
917 ScalarType const a13 = F[0] * F[1];
918 ScalarType const a14 = F[3] * F[4];
919 ScalarType const a15 = F[6] * F[7];
920 ScalarType const a16 = a13 * lambda + mu * (2 * a13 + a14 + a15);
921 ScalarType const a17 = F[0] * F[2];
922 ScalarType const a18 = F[3] * F[5];
923 ScalarType const a19 = F[6] * F[8];
924 ScalarType const a20 = a17 * lambda + mu * (2 * a17 + a18 + a19);
925 ScalarType const a21 = F[0] * F[3];
926 ScalarType const a22 = F[1] * F[4];
927 ScalarType const a23 = F[2] * F[5];
928 ScalarType const a24 = a21 * lambda + mu * (2 * a21 + a22 + a23);
929 ScalarType const a25 = lambda * F[0];
930 ScalarType const a26 = mu * F[3];
931 ScalarType const a27 = a25 * F[4] + a26 * F[1];
932 ScalarType const a28 = a25 * F[5] + a26 * F[2];
933 ScalarType const a29 = F[0] * F[6];
934 ScalarType const a30 = F[1] * F[7];
935 ScalarType const a31 = F[2] * F[8];
936 ScalarType const a32 = a29 * lambda + mu * (2 * a29 + a30 + a31);
937 ScalarType const a33 = mu * F[6];
938 ScalarType const a34 = a25 * F[7] + a33 * F[1];
939 ScalarType const a35 = a25 * F[8] + a33 * F[2];
940 ScalarType const a36 = a0 + a8;
941 ScalarType const a37 = F[1] * F[2];
942 ScalarType const a38 = F[4] * F[5];
943 ScalarType const a39 = F[7] * F[8];
944 ScalarType const a40 = a37 * lambda + mu * (2 * a37 + a38 + a39);
945 ScalarType const a41 = lambda * F[1];
946 ScalarType const a42 = mu * F[4];
947 ScalarType const a43 = a41 * F[3] + a42 * F[0];
948 ScalarType const a44 = a22 * lambda + mu * (a21 + 2 * a22 + a23);
949 ScalarType const a45 = a41 * F[5] + a42 * F[2];
950 ScalarType const a46 = mu * F[7];
951 ScalarType const a47 = a41 * F[6] + a46 * F[0];
952 ScalarType const a48 = a30 * lambda + mu * (a29 + 2 * a30 + a31);
953 ScalarType const a49 = a41 * F[8] + a46 * F[2];
954 ScalarType const a50 = a9 - 1;
955 ScalarType const a51 = a0 + a11;
956 ScalarType const a52 = lambda * F[2];
957 ScalarType const a53 = mu * F[5];
958 ScalarType const a54 = a52 * F[3] + a53 * F[0];
959 ScalarType const a55 = a52 * F[4] + a53 * F[1];
960 ScalarType const a56 = a23 * lambda + mu * (a21 + a22 + 2 * a23);
961 ScalarType const a57 = mu * F[8];
962 ScalarType const a58 = a52 * F[6] + a57 * F[0];
963 ScalarType const a59 = a52 * F[7] + a57 * F[1];
964 ScalarType const a60 = a31 * lambda + mu * (a29 + a30 + 2 * a31);
965 ScalarType const a61 = a14 * lambda + mu * (a13 + 2 * a14 + a15);
966 ScalarType const a62 = a18 * lambda + mu * (a17 + 2 * a18 + a19);
967 ScalarType const a63 = F[3] * F[6];
968 ScalarType const a64 = F[4] * F[7];
969 ScalarType const a65 = F[5] * F[8];
970 ScalarType const a66 = a63 * lambda + mu * (2 * a63 + a64 + a65);
971 ScalarType const a67 = lambda * F[3];
972 ScalarType const a68 = a33 * F[4] + a67 * F[7];
973 ScalarType const a69 = a33 * F[5] + a67 * F[8];
974 ScalarType const a70 = a38 * lambda + mu * (a37 + 2 * a38 + a39);
975 ScalarType const a71 = lambda * F[4];
976 ScalarType const a72 = a26 * F[7] + a71 * F[6];
977 ScalarType const a73 = a64 * lambda + mu * (a63 + 2 * a64 + a65);
978 ScalarType const a74 = a46 * F[5] + a71 * F[8];
979 ScalarType const a75 = a11 + a8;
980 ScalarType const a76 = lambda * F[5];
981 ScalarType const a77 = a26 * F[8] + a76 * F[6];
982 ScalarType const a78 = a42 * F[8] + a76 * F[7];
983 ScalarType const a79 = a65 * lambda + mu * (a63 + a64 + 2 * a65);
984 ScalarType const a80 = a15 * lambda + mu * (a13 + a14 + 2 * a15);
985 ScalarType const a81 = a19 * lambda + mu * (a17 + a18 + 2 * a19);
986 ScalarType const a82 = a39 * lambda + mu * (a37 + a38 + 2 * a39);
987 H[0] = a0 * lambda + a12 + mu * (3 * a0 + a3 + a7);
988 H[1] = a16;
989 H[2] = a20;
990 H[3] = a24;
991 H[4] = a27;
992 H[5] = a28;
993 H[6] = a32;
994 H[7] = a34;
995 H[8] = a35;
996 H[9] = a16;
997 H[10] = a1 * lambda + a12 + mu * (3 * a1 + a10 + a36 + a6);
998 H[11] = a40;
999 H[12] = a43;
1000 H[13] = a44;
1001 H[14] = a45;
1002 H[15] = a47;
1003 H[16] = a48;
1004 H[17] = a49;
1005 H[18] = a20;
1006 H[19] = a40;
1007 H[20] = a12 + a5 * lambda + mu * (a1 + 3 * a5 + a50 + a51);
1008 H[21] = a54;
1009 H[22] = a55;
1010 H[23] = a56;
1011 H[24] = a58;
1012 H[25] = a59;
1013 H[26] = a60;
1014 H[27] = a24;
1015 H[28] = a43;
1016 H[29] = a54;
1017 H[30] = a12 + a2 * lambda + mu * (3 * a2 + a36 + a4 + a50);
1018 H[31] = a61;
1019 H[32] = a62;
1020 H[33] = a66;
1021 H[34] = a68;
1022 H[35] = a69;
1023 H[36] = a27;
1024 H[37] = a44;
1025 H[38] = a55;
1026 H[39] = a61;
1027 H[40] = a12 + a8 * lambda + mu * (a10 + a3 + a50 + 3 * a8);
1028 H[41] = a70;
1029 H[42] = a72;
1030 H[43] = a73;
1031 H[44] = a74;
1032 H[45] = a28;
1033 H[46] = a45;
1034 H[47] = a56;
1035 H[48] = a62;
1036 H[49] = a70;
1037 H[50] = a12 + a9 * lambda + mu * (a2 + a6 + a75 + 3 * a9);
1038 H[51] = a77;
1039 H[52] = a78;
1040 H[53] = a79;
1041 H[54] = a32;
1042 H[55] = a47;
1043 H[56] = a58;
1044 H[57] = a66;
1045 H[58] = a72;
1046 H[59] = a77;
1047 H[60] = a12 + a4 * lambda + mu * (a10 + a2 + 3 * a4 + a51 - 1);
1048 H[61] = a80;
1049 H[62] = a81;
1050 H[63] = a34;
1051 H[64] = a48;
1052 H[65] = a59;
1053 H[66] = a68;
1054 H[67] = a73;
1055 H[68] = a78;
1056 H[69] = a80;
1057 H[70] = a10 * lambda + a12 + mu * (a1 + 3 * a10 + a4 + a75 - 1);
1058 H[71] = a82;
1059 H[72] = a35;
1060 H[73] = a49;
1061 H[74] = a60;
1062 H[75] = a69;
1063 H[76] = a74;
1064 H[77] = a79;
1065 H[78] = a81;
1066 H[79] = a82;
1067 H[80] = a11 * lambda + a12 + mu * (a10 + 3 * a11 + a7 + a9);
1068 return H;
1069}
1070
1071template <
1072 math::linalg::mini::CReadableVectorizedMatrix TMatrix,
1073 math::linalg::mini::CWriteableVectorizedMatrix TMatrixGF>
1074PBAT_HOST_DEVICE typename TMatrix::ScalarType SaintVenantKirchhoffEnergy<3>::evalWithGrad(
1075 [[maybe_unused]] TMatrix const& F,
1076 [[maybe_unused]] typename TMatrix::ScalarType mu,
1077 [[maybe_unused]] typename TMatrix::ScalarType lambda,
1078 TMatrixGF& gF) const
1079{
1080 static_assert(
1081 TMatrixGF::kRows == 9 and TMatrixGF::kCols == 1,
1082 "Grad w.r.t. F must have dimensions 9x1");
1083 using ScalarType = typename TMatrix::ScalarType;
1084 ScalarType psi;
1085 ScalarType const a0 = (1.0 / 2.0) * ((F[0]) * (F[0])) + (1.0 / 2.0) * ((F[1]) * (F[1])) +
1086 (1.0 / 2.0) * ((F[2]) * (F[2]));
1087 ScalarType const a1 = (1.0 / 2.0) * ((F[3]) * (F[3])) + (1.0 / 2.0) * ((F[4]) * (F[4])) +
1088 (1.0 / 2.0) * ((F[5]) * (F[5]));
1089 ScalarType const a2 = (1.0 / 2.0) * ((F[6]) * (F[6])) + (1.0 / 2.0) * ((F[7]) * (F[7])) +
1090 (1.0 / 2.0) * ((F[8]) * (F[8]));
1091 ScalarType const a3 = a0 + a1 + a2 - 3.0 / 2.0;
1092 ScalarType const a4 = a0 - 1.0 / 2.0;
1093 ScalarType const a5 = a1 - 1.0 / 2.0;
1094 ScalarType const a6 = a2 - 1.0 / 2.0;
1095 ScalarType const a7 = (1.0 / 2.0) * F[0];
1096 ScalarType const a8 = (1.0 / 2.0) * F[1];
1097 ScalarType const a9 = (1.0 / 2.0) * F[2];
1098 ScalarType const a10 = a7 * F[3] + a8 * F[4] + a9 * F[5];
1099 ScalarType const a11 = a7 * F[6] + a8 * F[7] + a9 * F[8];
1100 ScalarType const a12 =
1101 (1.0 / 2.0) * F[3] * F[6] + (1.0 / 2.0) * F[4] * F[7] + (1.0 / 2.0) * F[5] * F[8];
1102 ScalarType const a13 = a3 * lambda;
1103 ScalarType const a14 = 2 * a4;
1104 ScalarType const a15 = 2 * a10;
1105 ScalarType const a16 = 2 * a11;
1106 ScalarType const a17 = 2 * a5;
1107 ScalarType const a18 = 2 * a12;
1108 ScalarType const a19 = 2 * a6;
1109 psi = (1.0 / 2.0) * ((a3) * (a3)) * lambda +
1110 mu * (2 * ((a10) * (a10)) + 2 * ((a11) * (a11)) + 2 * ((a12) * (a12)) + ((a4) * (a4)) +
1111 ((a5) * (a5)) + ((a6) * (a6)));
1112 gF[0] = a13 * F[0] + mu * (a14 * F[0] + a15 * F[3] + a16 * F[6]);
1113 gF[1] = a13 * F[1] + mu * (a14 * F[1] + a15 * F[4] + a16 * F[7]);
1114 gF[2] = a13 * F[2] + mu * (a14 * F[2] + a15 * F[5] + a16 * F[8]);
1115 gF[3] = a13 * F[3] + mu * (a15 * F[0] + a17 * F[3] + a18 * F[6]);
1116 gF[4] = a13 * F[4] + mu * (a15 * F[1] + a17 * F[4] + a18 * F[7]);
1117 gF[5] = a13 * F[5] + mu * (a15 * F[2] + a17 * F[5] + a18 * F[8]);
1118 gF[6] = a13 * F[6] + mu * (a16 * F[0] + a18 * F[3] + a19 * F[6]);
1119 gF[7] = a13 * F[7] + mu * (a16 * F[1] + a18 * F[4] + a19 * F[7]);
1120 gF[8] = a13 * F[8] + mu * (a16 * F[2] + a18 * F[5] + a19 * F[8]);
1121 return psi;
1122}
1123
1124template <
1125 math::linalg::mini::CReadableVectorizedMatrix TMatrix,
1126 math::linalg::mini::CWriteableVectorizedMatrix TMatrixGF,
1127 math::linalg::mini::CWriteableVectorizedMatrix TMatrixHF>
1128PBAT_HOST_DEVICE typename TMatrix::ScalarType SaintVenantKirchhoffEnergy<3>::evalWithGradAndHessian(
1129 [[maybe_unused]] TMatrix const& F,
1130 [[maybe_unused]] typename TMatrix::ScalarType mu,
1131 [[maybe_unused]] typename TMatrix::ScalarType lambda,
1132 TMatrixGF& gF,
1133 TMatrixHF& HF) const
1134{
1135 static_assert(
1136 TMatrixGF::kRows == 9 and TMatrixGF::kCols == 1,
1137 "Grad w.r.t. F must have dimensions 9x1");
1138 static_assert(
1139 TMatrixHF::kRows == 9 and TMatrixHF::kCols == 9,
1140 "Hessian w.r.t. F must have dimensions 9x9");
1141 using ScalarType = typename TMatrix::ScalarType;
1142 ScalarType psi;
1143 ScalarType const a0 = ((F[0]) * (F[0]));
1144 ScalarType const a1 = ((F[1]) * (F[1]));
1145 ScalarType const a2 = ((F[2]) * (F[2]));
1146 ScalarType const a3 = (1.0 / 2.0) * a0 + (1.0 / 2.0) * a1 + (1.0 / 2.0) * a2;
1147 ScalarType const a4 = ((F[3]) * (F[3]));
1148 ScalarType const a5 = ((F[4]) * (F[4]));
1149 ScalarType const a6 = ((F[5]) * (F[5]));
1150 ScalarType const a7 = (1.0 / 2.0) * a4 + (1.0 / 2.0) * a5 + (1.0 / 2.0) * a6;
1151 ScalarType const a8 = ((F[6]) * (F[6]));
1152 ScalarType const a9 = ((F[7]) * (F[7]));
1153 ScalarType const a10 = ((F[8]) * (F[8]));
1154 ScalarType const a11 = (1.0 / 2.0) * a10 + (1.0 / 2.0) * a8 + (1.0 / 2.0) * a9;
1155 ScalarType const a12 = a11 + a3 + a7 - 3.0 / 2.0;
1156 ScalarType const a13 = a3 - 1.0 / 2.0;
1157 ScalarType const a14 = a7 - 1.0 / 2.0;
1158 ScalarType const a15 = a11 - 1.0 / 2.0;
1159 ScalarType const a16 = F[0] * F[3];
1160 ScalarType const a17 = F[1] * F[4];
1161 ScalarType const a18 = F[2] * F[5];
1162 ScalarType const a19 = (1.0 / 2.0) * a16 + (1.0 / 2.0) * a17 + (1.0 / 2.0) * a18;
1163 ScalarType const a20 = F[0] * F[6];
1164 ScalarType const a21 = F[1] * F[7];
1165 ScalarType const a22 = F[2] * F[8];
1166 ScalarType const a23 = (1.0 / 2.0) * a20 + (1.0 / 2.0) * a21 + (1.0 / 2.0) * a22;
1167 ScalarType const a24 = F[3] * F[6];
1168 ScalarType const a25 = F[4] * F[7];
1169 ScalarType const a26 = F[5] * F[8];
1170 ScalarType const a27 = (1.0 / 2.0) * a24 + (1.0 / 2.0) * a25 + (1.0 / 2.0) * a26;
1171 ScalarType const a28 = a12 * lambda;
1172 ScalarType const a29 = 2 * a13;
1173 ScalarType const a30 = 2 * a19;
1174 ScalarType const a31 = 2 * a23;
1175 ScalarType const a32 = 2 * a14;
1176 ScalarType const a33 = 2 * a27;
1177 ScalarType const a34 = 2 * a15;
1178 ScalarType const a35 = a1 + a4;
1179 ScalarType const a36 = a2 - 1;
1180 ScalarType const a37 = a36 + a8;
1181 ScalarType const a38 = F[0] * F[1];
1182 ScalarType const a39 = F[3] * F[4];
1183 ScalarType const a40 = F[6] * F[7];
1184 ScalarType const a41 = a38 * lambda + mu * (2 * a38 + a39 + a40);
1185 ScalarType const a42 = F[0] * F[2];
1186 ScalarType const a43 = F[3] * F[5];
1187 ScalarType const a44 = F[6] * F[8];
1188 ScalarType const a45 = a42 * lambda + mu * (2 * a42 + a43 + a44);
1189 ScalarType const a46 = a16 * lambda + mu * (2 * a16 + a17 + a18);
1190 ScalarType const a47 = lambda * F[0];
1191 ScalarType const a48 = mu * F[3];
1192 ScalarType const a49 = a47 * F[4] + a48 * F[1];
1193 ScalarType const a50 = a47 * F[5] + a48 * F[2];
1194 ScalarType const a51 = a20 * lambda + mu * (2 * a20 + a21 + a22);
1195 ScalarType const a52 = mu * F[6];
1196 ScalarType const a53 = a47 * F[7] + a52 * F[1];
1197 ScalarType const a54 = a47 * F[8] + a52 * F[2];
1198 ScalarType const a55 = a0 + a5;
1199 ScalarType const a56 = F[1] * F[2];
1200 ScalarType const a57 = F[4] * F[5];
1201 ScalarType const a58 = F[7] * F[8];
1202 ScalarType const a59 = a56 * lambda + mu * (2 * a56 + a57 + a58);
1203 ScalarType const a60 = lambda * F[1];
1204 ScalarType const a61 = mu * F[4];
1205 ScalarType const a62 = a60 * F[3] + a61 * F[0];
1206 ScalarType const a63 = a17 * lambda + mu * (a16 + 2 * a17 + a18);
1207 ScalarType const a64 = a60 * F[5] + a61 * F[2];
1208 ScalarType const a65 = mu * F[7];
1209 ScalarType const a66 = a60 * F[6] + a65 * F[0];
1210 ScalarType const a67 = a21 * lambda + mu * (a20 + 2 * a21 + a22);
1211 ScalarType const a68 = a60 * F[8] + a65 * F[2];
1212 ScalarType const a69 = a6 - 1;
1213 ScalarType const a70 = a0 + a10;
1214 ScalarType const a71 = lambda * F[2];
1215 ScalarType const a72 = mu * F[5];
1216 ScalarType const a73 = a71 * F[3] + a72 * F[0];
1217 ScalarType const a74 = a71 * F[4] + a72 * F[1];
1218 ScalarType const a75 = a18 * lambda + mu * (a16 + a17 + 2 * a18);
1219 ScalarType const a76 = mu * F[8];
1220 ScalarType const a77 = a71 * F[6] + a76 * F[0];
1221 ScalarType const a78 = a71 * F[7] + a76 * F[1];
1222 ScalarType const a79 = a22 * lambda + mu * (a20 + a21 + 2 * a22);
1223 ScalarType const a80 = a39 * lambda + mu * (a38 + 2 * a39 + a40);
1224 ScalarType const a81 = a43 * lambda + mu * (a42 + 2 * a43 + a44);
1225 ScalarType const a82 = a24 * lambda + mu * (2 * a24 + a25 + a26);
1226 ScalarType const a83 = lambda * F[3];
1227 ScalarType const a84 = a52 * F[4] + a83 * F[7];
1228 ScalarType const a85 = a52 * F[5] + a83 * F[8];
1229 ScalarType const a86 = a57 * lambda + mu * (a56 + 2 * a57 + a58);
1230 ScalarType const a87 = lambda * F[4];
1231 ScalarType const a88 = a48 * F[7] + a87 * F[6];
1232 ScalarType const a89 = a25 * lambda + mu * (a24 + 2 * a25 + a26);
1233 ScalarType const a90 = a65 * F[5] + a87 * F[8];
1234 ScalarType const a91 = a10 + a5;
1235 ScalarType const a92 = lambda * F[5];
1236 ScalarType const a93 = a48 * F[8] + a92 * F[6];
1237 ScalarType const a94 = a61 * F[8] + a92 * F[7];
1238 ScalarType const a95 = a26 * lambda + mu * (a24 + a25 + 2 * a26);
1239 ScalarType const a96 = a40 * lambda + mu * (a38 + a39 + 2 * a40);
1240 ScalarType const a97 = a44 * lambda + mu * (a42 + a43 + 2 * a44);
1241 ScalarType const a98 = a58 * lambda + mu * (a56 + a57 + 2 * a58);
1242 psi = (1.0 / 2.0) * ((a12) * (a12)) * lambda +
1243 mu * (((a13) * (a13)) + ((a14) * (a14)) + ((a15) * (a15)) + 2 * ((a19) * (a19)) +
1244 2 * ((a23) * (a23)) + 2 * ((a27) * (a27)));
1245 gF[0] = a28 * F[0] + mu * (a29 * F[0] + a30 * F[3] + a31 * F[6]);
1246 gF[1] = a28 * F[1] + mu * (a29 * F[1] + a30 * F[4] + a31 * F[7]);
1247 gF[2] = a28 * F[2] + mu * (a29 * F[2] + a30 * F[5] + a31 * F[8]);
1248 gF[3] = a28 * F[3] + mu * (a30 * F[0] + a32 * F[3] + a33 * F[6]);
1249 gF[4] = a28 * F[4] + mu * (a30 * F[1] + a32 * F[4] + a33 * F[7]);
1250 gF[5] = a28 * F[5] + mu * (a30 * F[2] + a32 * F[5] + a33 * F[8]);
1251 gF[6] = a28 * F[6] + mu * (a31 * F[0] + a33 * F[3] + a34 * F[6]);
1252 gF[7] = a28 * F[7] + mu * (a31 * F[1] + a33 * F[4] + a34 * F[7]);
1253 gF[8] = a28 * F[8] + mu * (a31 * F[2] + a33 * F[5] + a34 * F[8]);
1254 HF[0] = a0 * lambda + a28 + mu * (3 * a0 + a35 + a37);
1255 HF[1] = a41;
1256 HF[2] = a45;
1257 HF[3] = a46;
1258 HF[4] = a49;
1259 HF[5] = a50;
1260 HF[6] = a51;
1261 HF[7] = a53;
1262 HF[8] = a54;
1263 HF[9] = a41;
1264 HF[10] = a1 * lambda + a28 + mu * (3 * a1 + a36 + a55 + a9);
1265 HF[11] = a59;
1266 HF[12] = a62;
1267 HF[13] = a63;
1268 HF[14] = a64;
1269 HF[15] = a66;
1270 HF[16] = a67;
1271 HF[17] = a68;
1272 HF[18] = a45;
1273 HF[19] = a59;
1274 HF[20] = a2 * lambda + a28 + mu * (a1 + 3 * a2 + a69 + a70);
1275 HF[21] = a73;
1276 HF[22] = a74;
1277 HF[23] = a75;
1278 HF[24] = a77;
1279 HF[25] = a78;
1280 HF[26] = a79;
1281 HF[27] = a46;
1282 HF[28] = a62;
1283 HF[29] = a73;
1284 HF[30] = a28 + a4 * lambda + mu * (3 * a4 + a55 + a69 + a8);
1285 HF[31] = a80;
1286 HF[32] = a81;
1287 HF[33] = a82;
1288 HF[34] = a84;
1289 HF[35] = a85;
1290 HF[36] = a49;
1291 HF[37] = a63;
1292 HF[38] = a74;
1293 HF[39] = a80;
1294 HF[40] = a28 + a5 * lambda + mu * (a35 + 3 * a5 + a69 + a9);
1295 HF[41] = a86;
1296 HF[42] = a88;
1297 HF[43] = a89;
1298 HF[44] = a90;
1299 HF[45] = a50;
1300 HF[46] = a64;
1301 HF[47] = a75;
1302 HF[48] = a81;
1303 HF[49] = a86;
1304 HF[50] = a28 + a6 * lambda + mu * (a36 + a4 + 3 * a6 + a91);
1305 HF[51] = a93;
1306 HF[52] = a94;
1307 HF[53] = a95;
1308 HF[54] = a51;
1309 HF[55] = a66;
1310 HF[56] = a77;
1311 HF[57] = a82;
1312 HF[58] = a88;
1313 HF[59] = a93;
1314 HF[60] = a28 + a8 * lambda + mu * (a4 + a70 + 3 * a8 + a9 - 1);
1315 HF[61] = a96;
1316 HF[62] = a97;
1317 HF[63] = a53;
1318 HF[64] = a67;
1319 HF[65] = a78;
1320 HF[66] = a84;
1321 HF[67] = a89;
1322 HF[68] = a94;
1323 HF[69] = a96;
1324 HF[70] = a28 + a9 * lambda + mu * (a1 + a8 + 3 * a9 + a91 - 1);
1325 HF[71] = a98;
1326 HF[72] = a54;
1327 HF[73] = a68;
1328 HF[74] = a79;
1329 HF[75] = a85;
1330 HF[76] = a90;
1331 HF[77] = a95;
1332 HF[78] = a97;
1333 HF[79] = a98;
1334 HF[80] = a10 * lambda + a28 + mu * (3 * a10 + a37 + a6 + a9);
1335 return psi;
1336}
1337
1338template <
1339 math::linalg::mini::CReadableVectorizedMatrix TMatrix,
1340 math::linalg::mini::CWriteableVectorizedMatrix TMatrixGF,
1341 math::linalg::mini::CWriteableVectorizedMatrix TMatrixHF>
1342PBAT_HOST_DEVICE void SaintVenantKirchhoffEnergy<3>::gradAndHessian(
1343 [[maybe_unused]] TMatrix const& F,
1344 [[maybe_unused]] typename TMatrix::ScalarType mu,
1345 [[maybe_unused]] typename TMatrix::ScalarType lambda,
1346 TMatrixGF& gF,
1347 TMatrixHF& HF) const
1348{
1349 static_assert(
1350 TMatrixGF::kRows == 9 and TMatrixGF::kCols == 1,
1351 "Grad w.r.t. F must have dimensions 9x1");
1352 static_assert(
1353 TMatrixHF::kRows == 9 and TMatrixHF::kCols == 9,
1354 "Hessian w.r.t. F must have dimensions 9x9");
1355 using ScalarType = typename TMatrix::ScalarType;
1356 ScalarType const a0 = ((F[0]) * (F[0]));
1357 ScalarType const a1 = ((F[1]) * (F[1]));
1358 ScalarType const a2 = ((F[2]) * (F[2]));
1359 ScalarType const a3 = (1.0 / 2.0) * a0 + (1.0 / 2.0) * a1 + (1.0 / 2.0) * a2;
1360 ScalarType const a4 = ((F[3]) * (F[3]));
1361 ScalarType const a5 = ((F[4]) * (F[4]));
1362 ScalarType const a6 = ((F[5]) * (F[5]));
1363 ScalarType const a7 = (1.0 / 2.0) * a4 + (1.0 / 2.0) * a5 + (1.0 / 2.0) * a6;
1364 ScalarType const a8 = ((F[6]) * (F[6]));
1365 ScalarType const a9 = ((F[7]) * (F[7]));
1366 ScalarType const a10 = ((F[8]) * (F[8]));
1367 ScalarType const a11 = (1.0 / 2.0) * a10 + (1.0 / 2.0) * a8 + (1.0 / 2.0) * a9;
1368 ScalarType const a12 = lambda * (a11 + a3 + a7 - 3.0 / 2.0);
1369 ScalarType const a13 = 2 * a3 - 1;
1370 ScalarType const a14 = F[0] * F[3];
1371 ScalarType const a15 = F[1] * F[4];
1372 ScalarType const a16 = F[2] * F[5];
1373 ScalarType const a17 = a14 + a15 + a16;
1374 ScalarType const a18 = F[0] * F[6];
1375 ScalarType const a19 = F[1] * F[7];
1376 ScalarType const a20 = F[2] * F[8];
1377 ScalarType const a21 = a18 + a19 + a20;
1378 ScalarType const a22 = 2 * a7 - 1;
1379 ScalarType const a23 = F[3] * F[6];
1380 ScalarType const a24 = F[4] * F[7];
1381 ScalarType const a25 = F[5] * F[8];
1382 ScalarType const a26 = a23 + a24 + a25;
1383 ScalarType const a27 = 2 * a11 - 1;
1384 ScalarType const a28 = a1 + a4;
1385 ScalarType const a29 = a2 - 1;
1386 ScalarType const a30 = a29 + a8;
1387 ScalarType const a31 = F[0] * F[1];
1388 ScalarType const a32 = F[3] * F[4];
1389 ScalarType const a33 = F[6] * F[7];
1390 ScalarType const a34 = a31 * lambda + mu * (2 * a31 + a32 + a33);
1391 ScalarType const a35 = F[0] * F[2];
1392 ScalarType const a36 = F[3] * F[5];
1393 ScalarType const a37 = F[6] * F[8];
1394 ScalarType const a38 = a35 * lambda + mu * (2 * a35 + a36 + a37);
1395 ScalarType const a39 = a14 * lambda + mu * (2 * a14 + a15 + a16);
1396 ScalarType const a40 = lambda * F[0];
1397 ScalarType const a41 = mu * F[3];
1398 ScalarType const a42 = a40 * F[4] + a41 * F[1];
1399 ScalarType const a43 = a40 * F[5] + a41 * F[2];
1400 ScalarType const a44 = a18 * lambda + mu * (2 * a18 + a19 + a20);
1401 ScalarType const a45 = mu * F[6];
1402 ScalarType const a46 = a40 * F[7] + a45 * F[1];
1403 ScalarType const a47 = a40 * F[8] + a45 * F[2];
1404 ScalarType const a48 = a0 + a5;
1405 ScalarType const a49 = F[1] * F[2];
1406 ScalarType const a50 = F[4] * F[5];
1407 ScalarType const a51 = F[7] * F[8];
1408 ScalarType const a52 = a49 * lambda + mu * (2 * a49 + a50 + a51);
1409 ScalarType const a53 = lambda * F[1];
1410 ScalarType const a54 = mu * F[4];
1411 ScalarType const a55 = a53 * F[3] + a54 * F[0];
1412 ScalarType const a56 = a15 * lambda + mu * (a14 + 2 * a15 + a16);
1413 ScalarType const a57 = a53 * F[5] + a54 * F[2];
1414 ScalarType const a58 = mu * F[7];
1415 ScalarType const a59 = a53 * F[6] + a58 * F[0];
1416 ScalarType const a60 = a19 * lambda + mu * (a18 + 2 * a19 + a20);
1417 ScalarType const a61 = a53 * F[8] + a58 * F[2];
1418 ScalarType const a62 = a6 - 1;
1419 ScalarType const a63 = a0 + a10;
1420 ScalarType const a64 = lambda * F[2];
1421 ScalarType const a65 = mu * F[5];
1422 ScalarType const a66 = a64 * F[3] + a65 * F[0];
1423 ScalarType const a67 = a64 * F[4] + a65 * F[1];
1424 ScalarType const a68 = a16 * lambda + mu * (a14 + a15 + 2 * a16);
1425 ScalarType const a69 = mu * F[8];
1426 ScalarType const a70 = a64 * F[6] + a69 * F[0];
1427 ScalarType const a71 = a64 * F[7] + a69 * F[1];
1428 ScalarType const a72 = a20 * lambda + mu * (a18 + a19 + 2 * a20);
1429 ScalarType const a73 = a32 * lambda + mu * (a31 + 2 * a32 + a33);
1430 ScalarType const a74 = a36 * lambda + mu * (a35 + 2 * a36 + a37);
1431 ScalarType const a75 = a23 * lambda + mu * (2 * a23 + a24 + a25);
1432 ScalarType const a76 = lambda * F[3];
1433 ScalarType const a77 = a45 * F[4] + a76 * F[7];
1434 ScalarType const a78 = a45 * F[5] + a76 * F[8];
1435 ScalarType const a79 = a50 * lambda + mu * (a49 + 2 * a50 + a51);
1436 ScalarType const a80 = lambda * F[4];
1437 ScalarType const a81 = a41 * F[7] + a80 * F[6];
1438 ScalarType const a82 = a24 * lambda + mu * (a23 + 2 * a24 + a25);
1439 ScalarType const a83 = a58 * F[5] + a80 * F[8];
1440 ScalarType const a84 = a10 + a5;
1441 ScalarType const a85 = lambda * F[5];
1442 ScalarType const a86 = a41 * F[8] + a85 * F[6];
1443 ScalarType const a87 = a54 * F[8] + a85 * F[7];
1444 ScalarType const a88 = a25 * lambda + mu * (a23 + a24 + 2 * a25);
1445 ScalarType const a89 = a33 * lambda + mu * (a31 + a32 + 2 * a33);
1446 ScalarType const a90 = a37 * lambda + mu * (a35 + a36 + 2 * a37);
1447 ScalarType const a91 = a51 * lambda + mu * (a49 + a50 + 2 * a51);
1448 gF[0] = a12 * F[0] + mu * (a13 * F[0] + a17 * F[3] + a21 * F[6]);
1449 gF[1] = a12 * F[1] + mu * (a13 * F[1] + a17 * F[4] + a21 * F[7]);
1450 gF[2] = a12 * F[2] + mu * (a13 * F[2] + a17 * F[5] + a21 * F[8]);
1451 gF[3] = a12 * F[3] + mu * (a17 * F[0] + a22 * F[3] + a26 * F[6]);
1452 gF[4] = a12 * F[4] + mu * (a17 * F[1] + a22 * F[4] + a26 * F[7]);
1453 gF[5] = a12 * F[5] + mu * (a17 * F[2] + a22 * F[5] + a26 * F[8]);
1454 gF[6] = a12 * F[6] + mu * (a21 * F[0] + a26 * F[3] + a27 * F[6]);
1455 gF[7] = a12 * F[7] + mu * (a21 * F[1] + a26 * F[4] + a27 * F[7]);
1456 gF[8] = a12 * F[8] + mu * (a21 * F[2] + a26 * F[5] + a27 * F[8]);
1457 HF[0] = a0 * lambda + a12 + mu * (3 * a0 + a28 + a30);
1458 HF[1] = a34;
1459 HF[2] = a38;
1460 HF[3] = a39;
1461 HF[4] = a42;
1462 HF[5] = a43;
1463 HF[6] = a44;
1464 HF[7] = a46;
1465 HF[8] = a47;
1466 HF[9] = a34;
1467 HF[10] = a1 * lambda + a12 + mu * (3 * a1 + a29 + a48 + a9);
1468 HF[11] = a52;
1469 HF[12] = a55;
1470 HF[13] = a56;
1471 HF[14] = a57;
1472 HF[15] = a59;
1473 HF[16] = a60;
1474 HF[17] = a61;
1475 HF[18] = a38;
1476 HF[19] = a52;
1477 HF[20] = a12 + a2 * lambda + mu * (a1 + 3 * a2 + a62 + a63);
1478 HF[21] = a66;
1479 HF[22] = a67;
1480 HF[23] = a68;
1481 HF[24] = a70;
1482 HF[25] = a71;
1483 HF[26] = a72;
1484 HF[27] = a39;
1485 HF[28] = a55;
1486 HF[29] = a66;
1487 HF[30] = a12 + a4 * lambda + mu * (3 * a4 + a48 + a62 + a8);
1488 HF[31] = a73;
1489 HF[32] = a74;
1490 HF[33] = a75;
1491 HF[34] = a77;
1492 HF[35] = a78;
1493 HF[36] = a42;
1494 HF[37] = a56;
1495 HF[38] = a67;
1496 HF[39] = a73;
1497 HF[40] = a12 + a5 * lambda + mu * (a28 + 3 * a5 + a62 + a9);
1498 HF[41] = a79;
1499 HF[42] = a81;
1500 HF[43] = a82;
1501 HF[44] = a83;
1502 HF[45] = a43;
1503 HF[46] = a57;
1504 HF[47] = a68;
1505 HF[48] = a74;
1506 HF[49] = a79;
1507 HF[50] = a12 + a6 * lambda + mu * (a29 + a4 + 3 * a6 + a84);
1508 HF[51] = a86;
1509 HF[52] = a87;
1510 HF[53] = a88;
1511 HF[54] = a44;
1512 HF[55] = a59;
1513 HF[56] = a70;
1514 HF[57] = a75;
1515 HF[58] = a81;
1516 HF[59] = a86;
1517 HF[60] = a12 + a8 * lambda + mu * (a4 + a63 + 3 * a8 + a9 - 1);
1518 HF[61] = a89;
1519 HF[62] = a90;
1520 HF[63] = a46;
1521 HF[64] = a60;
1522 HF[65] = a71;
1523 HF[66] = a77;
1524 HF[67] = a82;
1525 HF[68] = a87;
1526 HF[69] = a89;
1527 HF[70] = a12 + a9 * lambda + mu * (a1 + a8 + a84 + 3 * a9 - 1);
1528 HF[71] = a91;
1529 HF[72] = a47;
1530 HF[73] = a61;
1531 HF[74] = a72;
1532 HF[75] = a78;
1533 HF[76] = a83;
1534 HF[77] = a88;
1535 HF[78] = a90;
1536 HF[79] = a91;
1537 HF[80] = a10 * lambda + a12 + mu * (3 * a10 + a30 + a6 + a9);
1538}
1539
1540} // namespace physics
1541} // namespace pbat
1542
1543#endif // PBAT_PHYSICS_SAINTVENANTKIRCHHOFFENERGY_H
pbat::math::linalg::mini::SMatrix
Definition Matrix.h:121
pbat
The main namespace of the library.
Definition Aliases.h:15
pbat::physics::SaintVenantKirchhoffEnergy< 1 >::kDims
static auto constexpr kDims
Dimension of the space.
Definition SaintVenantKirchhoffEnergy.h:40
pbat::physics::SaintVenantKirchhoffEnergy< 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 SaintVenantKirchhoffEnergy.h:213
pbat::physics::SaintVenantKirchhoffEnergy< 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 SaintVenantKirchhoffEnergy.h:267
pbat::physics::SaintVenantKirchhoffEnergy< 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 SaintVenantKirchhoffEnergy.h:236
pbat::physics::SaintVenantKirchhoffEnergy< 1 >::SVector
pbat::math::linalg::mini::SVector< TScalar, M > SVector
Scalar vector type.
Definition SaintVenantKirchhoffEnergy.h:38
pbat::physics::SaintVenantKirchhoffEnergy< 1 >::SMatrix
pbat::math::linalg::mini::SMatrix< TScalar, M, N > SMatrix
Scalar matrix type.
Definition SaintVenantKirchhoffEnergy.h:35
pbat::physics::SaintVenantKirchhoffEnergy< 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::SaintVenantKirchhoffEnergy< 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::SaintVenantKirchhoffEnergy< 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 SaintVenantKirchhoffEnergy.h:148
pbat::physics::SaintVenantKirchhoffEnergy< 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::SaintVenantKirchhoffEnergy< 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 SaintVenantKirchhoffEnergy.h:520
pbat::physics::SaintVenantKirchhoffEnergy< 2 >::SVector
pbat::math::linalg::mini::SVector< TScalar, M > SVector
Scalar vector type.
Definition SaintVenantKirchhoffEnergy.h:303
pbat::physics::SaintVenantKirchhoffEnergy< 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 SaintVenantKirchhoffEnergy.h:554
pbat::physics::SaintVenantKirchhoffEnergy< 2 >::SMatrix
pbat::math::linalg::mini::SMatrix< TScalar, M, N > SMatrix
Scalar matrix type.
Definition SaintVenantKirchhoffEnergy.h:300
pbat::physics::SaintVenantKirchhoffEnergy< 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::SaintVenantKirchhoffEnergy< 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 SaintVenantKirchhoffEnergy.h:413
pbat::physics::SaintVenantKirchhoffEnergy< 2 >::kDims
static auto constexpr kDims
Dimension of the space.
Definition SaintVenantKirchhoffEnergy.h:305
pbat::physics::SaintVenantKirchhoffEnergy< 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 SaintVenantKirchhoffEnergy.h:626
pbat::physics::SaintVenantKirchhoffEnergy< 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 SaintVenantKirchhoffEnergy.h:1074
pbat::physics::SaintVenantKirchhoffEnergy< 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 SaintVenantKirchhoffEnergy.h:1342
pbat::physics::SaintVenantKirchhoffEnergy< 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::SaintVenantKirchhoffEnergy< 3 >::kDims
static auto constexpr kDims
Dimension of the space.
Definition SaintVenantKirchhoffEnergy.h:701
pbat::physics::SaintVenantKirchhoffEnergy< 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 SaintVenantKirchhoffEnergy.h:1128
pbat::physics::SaintVenantKirchhoffEnergy< 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 SaintVenantKirchhoffEnergy.h:809
pbat::physics::SaintVenantKirchhoffEnergy< 3 >::SMatrix
pbat::math::linalg::mini::SMatrix< TScalar, M, N > SMatrix
Scalar matrix type.
Definition SaintVenantKirchhoffEnergy.h:696
pbat::physics::SaintVenantKirchhoffEnergy< 3 >::SVector
pbat::math::linalg::mini::SVector< TScalar, M > SVector
Scalar vector type.
Definition SaintVenantKirchhoffEnergy.h:699
pbat::physics::SaintVenantKirchhoffEnergy< 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::SaintVenantKirchhoffEnergy
Definition SaintVenantKirchhoffEnergy.h:23