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
HyperElasticPotential.h
Go to the documentation of this file.
1
10
11#ifndef PBAT_FEM_HYPERELASTICPOTENTIAL_H
12#define PBAT_FEM_HYPERELASTICPOTENTIAL_H
13
14#include "Concepts.h"
15#include "DeformationGradient.h"
16#include "PhysicsBasedAnimationToolkitExport.h"
17#include "pbat/Aliases.h"
18#include "pbat/common/Eigen.h"
19#include "pbat/math/linalg/FilterEigenvalues.h"
20#include "pbat/math/linalg/SparsityPattern.h"
21#include "pbat/math/linalg/mini/Eigen.h"
22#include "pbat/math/linalg/mini/Product.h"
23#include "pbat/physics/HyperElasticity.h"
24#include "pbat/profiling/Profiling.h"
25
26#include <Eigen/Eigenvalues>
27#include <exception>
28#include <fmt/core.h>
29#include <span>
30#include <string>
31#include <tbb/parallel_for.h>
32
33namespace pbat::fem {
34
38enum class EHyperElasticSpdCorrection : std::uint32_t {
39 None, // No projection
40 Projection, // Project element hessian to nearest (in the 2-norm) SPD matrix
41 Absolute, // Flip negative eigenvalue signs
42};
43
49PBAT_API math::linalg::EEigenvalueFilter ToEigenvalueFilter(EHyperElasticSpdCorrection mode);
50
54enum EElementElasticityComputationFlags : int {
55 Potential = 1 << 0, // Compute element elastic potential at quadrature points
56 Gradient = 1 << 1, // Compute element gradient vectors at quadrature points
57 Hessian = 1 << 2, // Compute element hessian matrices at quadrature points
58};
59
90template <
91 CElement TElement,
92 int Dims,
93 physics::CHyperElasticEnergy THyperElasticEnergy,
94 class TDerivedE,
95 class TDerivedeg,
96 class TDerivedwg,
97 class TDerivedGNeg,
98 class TDerivedmug,
99 class TDerivedlambdag,
100 class TDerivedx,
101 class TDerivedUg,
102 class TDerivedGg,
103 class TDerivedHg>
104void ToElementElasticity(
105 Eigen::DenseBase<TDerivedE> const& E,
106 Eigen::Index nNodes,
107 Eigen::DenseBase<TDerivedeg> const& eg,
108 Eigen::DenseBase<TDerivedwg> const& wg,
109 Eigen::MatrixBase<TDerivedGNeg> const& GNeg,
110 Eigen::DenseBase<TDerivedmug> const& mug,
111 Eigen::DenseBase<TDerivedlambdag> const& lambdag,
112 Eigen::MatrixBase<TDerivedx> const& x,
113 Eigen::PlainObjectBase<TDerivedUg>& Ug,
114 Eigen::PlainObjectBase<TDerivedGg>& Gg,
115 Eigen::PlainObjectBase<TDerivedHg>& Hg,
116 int eFlags = EElementElasticityComputationFlags::Potential,
117 EHyperElasticSpdCorrection eSpdCorrection = EHyperElasticSpdCorrection::None);
118
147template <
148 physics::CHyperElasticEnergy THyperElasticEnergy,
149 CMesh TMesh,
150 class TDerivedeg,
151 class TDerivedwg,
152 class TDerivedGNeg,
153 class TDerivedmug,
154 class TDerivedlambdag,
155 class TDerivedx,
156 class TDerivedUg,
157 class TDerivedGg,
158 class TDerivedHg>
159inline void ToElementElasticity(
160 TMesh const& mesh,
161 Eigen::DenseBase<TDerivedeg> const& eg,
162 Eigen::DenseBase<TDerivedwg> const& wg,
163 Eigen::MatrixBase<TDerivedGNeg> const& GNeg,
164 Eigen::DenseBase<TDerivedmug> const& mug,
165 Eigen::DenseBase<TDerivedlambdag> const& lambdag,
166 Eigen::MatrixBase<TDerivedx> const& x,
167 Eigen::PlainObjectBase<TDerivedUg>& Ug,
168 Eigen::PlainObjectBase<TDerivedGg>& Gg,
169 Eigen::PlainObjectBase<TDerivedHg>& Hg,
170 int eFlags = EElementElasticityComputationFlags::Potential,
171 EHyperElasticSpdCorrection eSpdCorrection = EHyperElasticSpdCorrection::None)
172{
173 ToElementElasticity<typename TMesh::ElementType, TMesh::kDims, THyperElasticEnergy>(
174 mesh.E,
175 static_cast<typename TMesh::IndexType>(mesh.X.cols()),
176 eg.derived(),
177 wg.derived(),
178 GNeg.derived(),
179 mug.derived(),
180 lambdag.derived(),
181 x.derived(),
182 Ug,
183 Gg,
184 Hg,
185 eFlags,
186 eSpdCorrection);
187}
188
206template <
207 CElement TElement,
208 int Dims,
209 class TDerivedE,
210 class TDerivedeg,
211 class TDerivedHg,
212 class TDerivedIn,
213 class TDerivedOut>
214void GemmHyperElastic(
215 Eigen::DenseBase<TDerivedE> const& E,
216 Eigen::Index nNodes,
217 Eigen::DenseBase<TDerivedeg> const& eg,
218 Eigen::MatrixBase<TDerivedHg> const& Hg,
219 Eigen::MatrixBase<TDerivedIn> const& X,
220 Eigen::DenseBase<TDerivedOut>& Y);
221
236template <CMesh TMesh, class TDerivedeg, class TDerivedHg, class TDerivedIn, class TDerivedOut>
237inline void GemmHyperElastic(
238 TMesh const& mesh,
239 Eigen::DenseBase<TDerivedeg> const& eg,
240 Eigen::MatrixBase<TDerivedHg> const& Hg,
241 Eigen::MatrixBase<TDerivedIn> const& X,
242 Eigen::DenseBase<TDerivedOut>& Y)
243{
244 GemmHyperElastic<typename TMesh::ElementType, TMesh::kDims>(
245 mesh.E,
246 static_cast<typename TMesh::IndexType>(mesh.X.cols()),
247 eg.derived(),
248 Hg.derived(),
249 X.derived(),
250 Y.derived());
251}
252
265template <
266 CElement TElement,
267 int Dims,
268 Eigen::StorageOptions Options,
269 class TDerivedE,
270 class TDerivedeg>
271auto ElasticHessianSparsity(
272 Eigen::DenseBase<TDerivedE> const& E,
273 Eigen::Index nNodes,
274 Eigen::DenseBase<TDerivedeg> const& eg)
275 -> math::linalg::SparsityPattern<typename TDerivedE::Scalar, Options>;
276
282template <Eigen::StorageOptions Options, CMesh TMesh>
283auto ElasticHessianSparsity(TMesh const& mesh)
284 -> math::linalg::SparsityPattern<typename TMesh::IndexType, Options>
285{
286 return ElasticHessianSparsity<typename TMesh::ElementType, TMesh::kDims, Options>(
287 mesh.E,
288 static_cast<typename TMesh::IndexType>(mesh.X.cols()),
289 mesh.eg.derived());
290}
291
307template <
308 CElement TElement,
309 int Dims,
310 Eigen::StorageOptions Options,
311 class TDerivedE,
312 class TDerivedeg,
313 class TDerivedHg>
314auto HyperElasticHessian(
315 Eigen::DenseBase<TDerivedE> const& E,
316 Eigen::Index nNodes,
317 Eigen::DenseBase<TDerivedeg> const& eg,
318 Eigen::MatrixBase<TDerivedHg> const& Hg)
319 -> Eigen::SparseMatrix<typename TDerivedHg::Scalar, Options, typename TDerivedE::Scalar>;
320
332template <Eigen::StorageOptions Options, CMesh TMesh, class TDerivedeg, class TDerivedHg>
333auto HyperElasticHessian(
334 TMesh const& mesh,
335 Eigen::DenseBase<TDerivedeg> const& eg,
336 Eigen::MatrixBase<TDerivedHg> const& Hg)
337 -> Eigen::SparseMatrix<typename TDerivedHg::Scalar, Options, typename TMesh::IndexType>
338{
339 return HyperElasticHessian<typename TMesh::ElementType, TMesh::kDims, Options>(
340 mesh.E,
341 static_cast<typename TMesh::IndexType>(mesh.X.cols()),
342 eg.derived(),
343 Hg.derived());
344}
345
361template <
362 CElement TElement,
363 int Dims,
364 class TDerivedE,
365 class TDerivedHg,
366 Eigen::StorageOptions Options,
367 class TDerivedH>
368void ToHyperElasticHessian(
369 Eigen::DenseBase<TDerivedE> const& E,
370 Eigen::Index nNodes,
371 Eigen::DenseBase<TDerivedHg> const& Hg,
372 math::linalg::SparsityPattern<typename TDerivedE::Scalar, Options> const& sparsity,
373 Eigen::SparseCompressedBase<TDerivedH>& H);
374
387template <CMesh TMesh, class TDerivedHg, Eigen::StorageOptions Options, class TDerivedH>
388void ToHyperElasticHessian(
389 TMesh const& mesh,
390 Eigen::DenseBase<TDerivedHg> const& Hg,
391 math::linalg::SparsityPattern<typename TMesh::IndexType, Options> const& sparsity,
392 Eigen::SparseCompressedBase<typename TDerivedHg::Scalar>& H)
393{
394 ToHyperElasticHessian<typename TMesh::ElementType, TMesh::kDims, TDerivedHg, Options>(
395 mesh.E,
396 static_cast<typename TMesh::IndexType>(mesh.X.cols()),
397 Hg.derived(),
398 sparsity,
399 H.derived());
400}
401
402template <
403 CElement TElement,
404 int Dims,
405 Eigen::StorageOptions Options,
406 class TDerivedE,
407 class TDerivedHg>
408auto HyperElasticHessian(
409 Eigen::DenseBase<TDerivedE> const& E,
410 Eigen::Index nNodes,
411 Eigen::DenseBase<TDerivedHg> const& Hg,
412 math::linalg::SparsityPattern<typename TDerivedE::Scalar, Options> const& sparsity)
413 -> Eigen::SparseMatrix<typename TDerivedHg::Scalar, Options, typename TDerivedE::Scalar>;
414
425template <Eigen::StorageOptions Options, CMesh TMesh, class TDerivedHg>
426auto HyperElasticHessian(
427 TMesh const& mesh,
428 Eigen::DenseBase<TDerivedHg> const& Hg,
429 math::linalg::SparsityPattern<typename TMesh::IndexType, Options> const& sparsity)
430 -> Eigen::SparseMatrix<typename TDerivedHg::Scalar, Options, typename TMesh::IndexType>
431{
432 return HyperElasticHessian<typename TMesh::ElementType, TMesh::kDims, Options>(
433 mesh.E,
434 static_cast<typename TMesh::IndexType>(mesh.X.cols()),
435 Hg.derived(),
436 sparsity);
437}
438
454template <
455 CElement TElement,
456 int Dims,
457 class TDerivedE,
458 class TDerivedeg,
459 class TDerivedGg,
460 class TDerivedG>
461void ToHyperElasticGradient(
462 Eigen::DenseBase<TDerivedE> const& E,
463 Eigen::Index nNodes,
464 Eigen::DenseBase<TDerivedeg> const& eg,
465 Eigen::MatrixBase<TDerivedGg> const& Gg,
466 Eigen::PlainObjectBase<TDerivedG>& G);
467
494template <
495 CElement TElement,
496 int Dims,
497 physics::CHyperElasticEnergy THyperElasticEnergy,
498 class TDerivedE,
499 class TDerivedeg,
500 class TDerivedwg,
501 class TDerivedGNeg,
502 class TDerivedmug,
503 class TDerivedlambdag,
504 class TDerivedx,
505 class TDerivedGg,
506 class TDerivedG>
507inline void ToHyperElasticGradient(
508 Eigen::DenseBase<TDerivedE> const& E,
509 Eigen::Index nNodes,
510 Eigen::DenseBase<TDerivedeg> const& eg,
511 Eigen::DenseBase<TDerivedwg> const& wg,
512 Eigen::MatrixBase<TDerivedGNeg> const& GNeg,
513 Eigen::DenseBase<TDerivedmug> const& mug,
514 Eigen::DenseBase<TDerivedlambdag> const& lambdag,
515 Eigen::MatrixBase<TDerivedx> const& x,
516 Eigen::MatrixBase<TDerivedGg> const& Gg,
517 Eigen::PlainObjectBase<TDerivedG>& G)
518{
519 using ScalarType = typename TDerivedx::Scalar;
520 Eigen::Matrix<ScalarType, 0, 0> dummyUg, dummyHg;
521 ToElementElasticity<TElement, Dims, THyperElasticEnergy>(
522 E.derived(),
523 nNodes,
524 eg.derived(),
525 wg.derived(),
526 GNeg.derived(),
527 mug.derived(),
528 lambdag.derived(),
529 x.derived(),
530 dummyUg,
531 Gg.derived(),
532 dummyHg,
533 EElementElasticityComputationFlags::Gradient,
534 EHyperElasticSpdCorrection::None);
535 ToHyperElasticGradient<TElement, Dims>(
536 E.derived(),
537 nNodes,
538 eg.derived(),
539 Gg.derived(),
540 G.derived());
541}
542
555template <CMesh TMesh, class TDerivedeg, class TDerivedGg, class TDerivedOut>
556inline void ToHyperElasticGradient(
557 TMesh const& mesh,
558 Eigen::DenseBase<TDerivedeg> const& eg,
559 Eigen::MatrixBase<TDerivedGg> const& Gg,
560 Eigen::PlainObjectBase<TDerivedOut>& G)
561{
562 ToHyperElasticGradient<typename TMesh::ElementType, TMesh::kDims>(
563 mesh.E,
564 static_cast<typename TMesh::IndexType>(mesh.X.cols()),
565 eg.derived(),
566 Gg.derived(),
567 G.derived());
568}
569
595template <
596 CMesh TMesh,
597 physics::CHyperElasticEnergy THyperElasticEnergy,
598 class TDerivedE,
599 class TDerivedeg,
600 class TDerivedwg,
601 class TDerivedGNeg,
602 class TDerivedmug,
603 class TDerivedlambdag,
604 class TDerivedx,
605 class TDerivedGg,
606 class TDerivedG>
607inline void ToHyperElasticGradient(
608 TMesh const& mesh,
609 Eigen::DenseBase<TDerivedeg> const& eg,
610 Eigen::DenseBase<TDerivedwg> const& wg,
611 Eigen::MatrixBase<TDerivedGNeg> const& GNeg,
612 Eigen::DenseBase<TDerivedmug> const& mug,
613 Eigen::DenseBase<TDerivedlambdag> const& lambdag,
614 Eigen::MatrixBase<TDerivedx> const& x,
615 Eigen::MatrixBase<TDerivedGg> const& Gg,
616 Eigen::PlainObjectBase<TDerivedG>& G)
617{
618 using ScalarType = typename TDerivedx::Scalar;
619 Eigen::Matrix<ScalarType, 0, 0> dummyUg, dummyHg;
620 ToElementElasticity<THyperElasticEnergy>(
621 mesh,
622 eg.derived(),
623 wg.derived(),
624 GNeg.derived(),
625 mug.derived(),
626 lambdag.derived(),
627 x.derived(),
628 dummyUg,
629 Gg.derived(),
630 dummyHg,
631 EElementElasticityComputationFlags::Gradient,
632 EHyperElasticSpdCorrection::None);
633 ToHyperElasticGradient(mesh, eg.derived(), Gg.derived(), G.derived());
634}
635
650template <CElement TElement, int Dims, class TDerivedE, class TDerivedeg, class TDerivedGg>
651auto HyperElasticGradient(
652 Eigen::DenseBase<TDerivedE> const& E,
653 Eigen::Index nNodes,
654 Eigen::DenseBase<TDerivedeg> const& eg,
655 Eigen::MatrixBase<TDerivedGg> const& Gg)
656 -> Eigen::Vector<typename TDerivedGg::Scalar, Eigen::Dynamic>
657{
658 using ScalarType = typename TDerivedGg::Scalar;
659 auto const numberOfDofs = Dims * nNodes;
660 Eigen::Vector<ScalarType, Eigen::Dynamic> G(numberOfDofs);
661 ToHyperElasticGradient<TElement, Dims>(
662 E.derived(),
663 nNodes,
664 eg.derived(),
665 Gg.derived(),
666 G.derived());
667 return G;
668}
669
682template <CMesh TMesh, class TDerivedeg, class TDerivedGg>
683auto HyperElasticGradient(
684 TMesh const& mesh,
685 Eigen::DenseBase<TDerivedeg> const& eg,
686 Eigen::MatrixBase<TDerivedGg> const& Gg)
687{
688 return HyperElasticGradient<typename TMesh::ElementType, TMesh::kDims>(
689 mesh.E,
690 static_cast<typename TMesh::IndexType>(mesh.X.cols()),
691 eg.derived(),
692 Gg.derived());
693}
694
720template <
721 CElement TElement,
722 int Dims,
723 physics::CHyperElasticEnergy THyperElasticEnergy,
724 class TDerivedE,
725 class TDerivedeg,
726 class TDerivedwg,
727 class TDerivedGNeg,
728 class TDerivedmug,
729 class TDerivedlambdag,
730 class TDerivedx,
731 class TDerivedGg>
732auto HyperElasticGradient(
733 Eigen::DenseBase<TDerivedE> const& E,
734 Eigen::Index nNodes,
735 Eigen::DenseBase<TDerivedeg> const& eg,
736 Eigen::DenseBase<TDerivedwg> const& wg,
737 Eigen::MatrixBase<TDerivedGNeg> const& GNeg,
738 Eigen::DenseBase<TDerivedmug> const& mug,
739 Eigen::DenseBase<TDerivedlambdag> const& lambdag,
740 Eigen::MatrixBase<TDerivedx> const& x,
741 Eigen::MatrixBase<TDerivedGg> const& Gg)
742{
743 using ScalarType = typename TDerivedx::Scalar;
744 Eigen::Matrix<ScalarType, 0, 0> dummyUg, dummyHg;
745 ToElementElasticity<TElement, Dims, THyperElasticEnergy>(
746 E,
747 nNodes,
748 eg,
749 wg,
750 GNeg,
751 mug,
752 lambdag,
753 x,
754 dummyUg,
755 Gg,
756 dummyHg,
757 EElementElasticityComputationFlags::Gradient,
758 EHyperElasticSpdCorrection::None);
759 return HyperElasticGradient<TElement, Dims, THyperElasticEnergy>(E, nNodes, eg, Gg);
760}
761
784template <
785 CMesh TMesh,
786 physics::CHyperElasticEnergy THyperElasticEnergy,
787 class TDerivedeg,
788 class TDerivedwg,
789 class TDerivedGNeg,
790 class TDerivedmug,
791 class TDerivedlambdag,
792 class TDerivedx,
793 class TDerivedGg>
794auto HyperElasticGradient(
795 TMesh const& mesh,
796 Eigen::DenseBase<TDerivedeg> const& eg,
797 Eigen::DenseBase<TDerivedwg> const& wg,
798 Eigen::MatrixBase<TDerivedGNeg> const& GNeg,
799 Eigen::DenseBase<TDerivedmug> const& mug,
800 Eigen::DenseBase<TDerivedlambdag> const& lambdag,
801 Eigen::MatrixBase<TDerivedx> const& x,
802 Eigen::MatrixBase<TDerivedGg> const& Gg)
803{
804 return HyperElasticGradient<typename TMesh::ElementType, TMesh::kDims, THyperElasticEnergy>(
805 mesh.E,
806 static_cast<typename TMesh::IndexType>(mesh.X.cols()),
807 eg.derived(),
808 wg.derived(),
809 GNeg.derived(),
810 mug.derived(),
811 lambdag.derived(),
812 x.derived(),
813 Gg.derived());
814}
815
828template <
829 physics::CHyperElasticEnergy THyperElasticEnergy,
830 CMesh TMesh,
831 class TDerivedeg,
832 class TDerivedGg>
833auto HyperElasticGradient(
834 TMesh const& mesh,
835 Eigen::DenseBase<TDerivedeg> const& eg,
836 Eigen::MatrixBase<TDerivedGg> const& Gg)
837{
838 return HyperElasticGradient<typename TMesh::ElementType, TMesh::kDims, THyperElasticEnergy>(
839 mesh.E,
840 static_cast<typename TMesh::IndexType>(mesh.X.cols()),
841 eg.derived(),
842 Gg.derived());
843}
844
851template <class TDerivedUg>
852auto HyperElasticPotential(Eigen::DenseBase<TDerivedUg> const& Ug) -> typename TDerivedUg::Scalar
853{
854 return Ug.sum();
855}
856
879template <
880 CElement TElement,
881 int Dims,
882 physics::CHyperElasticEnergy THyperElasticEnergy,
883 class TDerivedE,
884 class TDerivedeg,
885 class TDerivedwg,
886 class TDerivedGNeg,
887 class TDerivedmug,
888 class TDerivedlambdag,
889 class TDerivedx>
890auto HyperElasticPotential(
891 Eigen::DenseBase<TDerivedE> const& E,
892 Eigen::Index nNodes,
893 Eigen::DenseBase<TDerivedeg> const& eg,
894 Eigen::DenseBase<TDerivedwg> const& wg,
895 Eigen::MatrixBase<TDerivedGNeg> const& GNeg,
896 Eigen::DenseBase<TDerivedmug> const& mug,
897 Eigen::DenseBase<TDerivedlambdag> const& lambdag,
898 Eigen::MatrixBase<TDerivedx> const& x)
899{
900 using ScalarType = typename TDerivedx::Scalar;
901 Eigen::Vector<ScalarType, Eigen::Dynamic> Ug;
902 Eigen::Matrix<ScalarType, 0, 0> dummyGg, dummyHg;
903 ToElementElasticity<TElement, Dims, THyperElasticEnergy>(
904 E,
905 nNodes,
906 eg,
907 wg,
908 GNeg,
909 mug,
910 lambdag,
911 x,
912 Ug,
913 dummyGg,
914 dummyHg,
915 EElementElasticityComputationFlags::Potential,
916 EHyperElasticSpdCorrection::None);
917 return HyperElasticPotential(Ug);
918}
919
920template <
921 CElement TElement,
922 int Dims,
923 physics::CHyperElasticEnergy THyperElasticEnergy,
924 class TDerivedE,
925 class TDerivedeg,
926 class TDerivedwg,
927 class TDerivedGNeg,
928 class TDerivedmug,
929 class TDerivedlambdag,
930 class TDerivedx,
931 class TDerivedUg,
932 class TDerivedGg,
933 class TDerivedHg>
934void ToElementElasticity(
935 Eigen::DenseBase<TDerivedE> const& E,
936 Eigen::Index nNodes,
937 Eigen::DenseBase<TDerivedeg> const& eg,
938 Eigen::DenseBase<TDerivedwg> const& wg,
939 Eigen::MatrixBase<TDerivedGNeg> const& GNeg,
940 Eigen::DenseBase<TDerivedmug> const& mug,
941 Eigen::DenseBase<TDerivedlambdag> const& lambdag,
942 Eigen::MatrixBase<TDerivedx> const& x,
943 Eigen::PlainObjectBase<TDerivedUg>& Ug,
944 Eigen::PlainObjectBase<TDerivedGg>& Gg,
945 Eigen::PlainObjectBase<TDerivedHg>& Hg,
946 int eFlags,
947 EHyperElasticSpdCorrection eSpdCorrection)
948{
949 PBAT_PROFILE_NAMED_SCOPE("pbat.fem.ToElementElasticity");
950
951 using ScalarType = typename TDerivedx::Scalar;
952
953 // Check inputs
954 if (x.size() != nNodes * Dims)
955 {
956 std::string const what = fmt::format(
957 "Generalized coordinate vector must have dimensions |#nodes|*kDims={}, but got "
958 "x.size()={}",
959 nNodes * Dims,
960 x.size());
961 throw std::invalid_argument(what);
962 }
963 if (GNeg.rows() != TElement::kNodes or GNeg.cols() != wg.size() * Dims)
964 {
965 std::string const what = fmt::format(
966 "Shape function gradient matrix must have dimensions |#nodes per element| x "
967 "|kDims*#quad pts|={}x{}, but got {}x{}",
968 TElement::kNodes,
969 wg.size() * Dims,
970 GNeg.rows(),
971 GNeg.cols());
972 throw std::invalid_argument(what);
973 }
974 if (eg.size() != wg.size() or mug.size() != wg.size() or lambdag.size() != wg.size())
975 {
976 std::string const what = fmt::format(
977 "Element index, quadrature weight, and Lame coefficient vectors must have the same "
978 "size of |#quad pts|={}, but got eg.size()={}, mug.size()={}, lambdag.size()={}",
979 wg.size(),
980 eg.size(),
981 mug.size(),
982 lambdag.size());
983 throw std::invalid_argument(what);
984 }
985
986 auto const nQuadPts = wg.size();
987 using SizeType = std::remove_cvref_t<decltype(nQuadPts)>;
988 auto constexpr kNodesPerElement = TElement::kNodes;
989 auto constexpr kDofsPerElement = kNodesPerElement * Dims;
990 namespace mini = math::linalg::mini;
991 using mini::FromEigen;
992 using mini::ToEigen;
993
994 if (eFlags & EElementElasticityComputationFlags::Potential)
995 {
996 Ug.resize(nQuadPts);
997 Ug.setZero();
998 }
999 if (eFlags & EElementElasticityComputationFlags::Gradient)
1000 {
1001 Gg.resize(kDofsPerElement, nQuadPts);
1002 Gg.setZero();
1003 }
1004 if (eFlags & EElementElasticityComputationFlags::Hessian)
1005 {
1006 Hg.resize(kDofsPerElement, kDofsPerElement * nQuadPts);
1007 Hg.setZero();
1008 }
1009
1010 THyperElasticEnergy Psi{};
1011 if (eFlags == EElementElasticityComputationFlags::Potential)
1012 {
1013 tbb::parallel_for(SizeType{0}, nQuadPts, [&](auto g) {
1014 auto const e = eg(g);
1015 auto const nodes = E.col(e);
1016 auto const xe = x.reshaped(Dims, nNodes)(Eigen::placeholders::all, nodes);
1017 auto const GPeg = GNeg.template block<kNodesPerElement, Dims>(0, g * Dims);
1018 Eigen::Matrix<ScalarType, Dims, Dims> const F = xe * GPeg;
1019 auto vecF = FromEigen(F);
1020 auto psiF = Psi.eval(vecF, mug(g), lambdag(g));
1021 Ug(g) += wg(g) * psiF;
1022 });
1023 }
1024 else if (
1025 eFlags == (EElementElasticityComputationFlags::Potential |
1026 EElementElasticityComputationFlags::Gradient))
1027 {
1028 tbb::parallel_for(SizeType{0}, nQuadPts, [&](auto g) {
1029 auto const e = eg(g);
1030 auto const nodes = E.col(e);
1031 auto const xe = x.reshaped(Dims, nNodes)(Eigen::placeholders::all, nodes);
1032 auto const GPeg = GNeg.template block<kNodesPerElement, Dims>(0, g * Dims);
1033 Eigen::Matrix<ScalarType, Dims, Dims> const F = xe * GPeg;
1034 auto vecF = FromEigen(F);
1035 mini::SVector<ScalarType, Dims * Dims> gradPsiF;
1036 auto const psiF = Psi.evalWithGrad(vecF, mug(g), lambdag(g), gradPsiF);
1037 auto const GP = FromEigen(GPeg);
1038 auto const GPsix = GradientWrtDofs<TElement, Dims>(gradPsiF, GP);
1039 Ug(g) += wg(g) * psiF;
1040 Gg.col(g) += wg(g) * ToEigen(GPsix);
1041 });
1042 }
1043 else if (
1044 eFlags == (EElementElasticityComputationFlags::Potential |
1045 EElementElasticityComputationFlags::Hessian))
1046 {
1047 tbb::parallel_for(SizeType{0}, nQuadPts, [&](auto g) {
1048 auto const e = eg(g);
1049 auto const nodes = E.col(e);
1050 auto const xe = x.reshaped(Dims, nNodes)(Eigen::placeholders::all, nodes);
1051 auto const gradPhi = GNeg.template block<kNodesPerElement, Dims>(0, g * Dims);
1052 Eigen::Matrix<ScalarType, Dims, Dims> const F = xe * gradPhi;
1053 auto vecF = FromEigen(F);
1054 auto psiF = Psi.eval(vecF, mug(g), lambdag(g));
1055 auto const hessPsiF = Psi.hessian(vecF, mug(g), lambdag(g));
1056 Eigen::Matrix<ScalarType, Dims * Dims, Dims * Dims> hessPsiFCorrected;
1057 math::linalg::FilterEigenvalues(
1058 ToEigen(hessPsiF),
1059 ToEigenvalueFilter(eSpdCorrection),
1060 hessPsiFCorrected);
1061 auto const GP = FromEigen(gradPhi);
1062 auto HPsix = HessianWrtDofs<TElement, Dims>(FromEigen(hessPsiFCorrected), GP);
1063 auto heg = Hg.template block<kDofsPerElement, kDofsPerElement>(0, g * kDofsPerElement);
1064 Ug(g) += wg(g) * psiF;
1065 heg += wg(g) * ToEigen(HPsix);
1066 });
1067 }
1068 else if (
1069 eFlags == (EElementElasticityComputationFlags::Potential |
1070 EElementElasticityComputationFlags::Gradient |
1071 EElementElasticityComputationFlags::Hessian))
1072 {
1073 tbb::parallel_for(SizeType{0}, nQuadPts, [&](auto g) {
1074 auto const e = eg(g);
1075 auto const nodes = E.col(e);
1076 Eigen::Matrix<ScalarType, Dims, kNodesPerElement> const xe =
1077 x.reshaped(Dims, nNodes)(Eigen::placeholders::all, nodes);
1078 auto const GPeg = GNeg.template block<kNodesPerElement, Dims>(0, g * Dims);
1079 Eigen::Matrix<ScalarType, Dims, Dims> const F = xe * GPeg;
1080 auto vecF = FromEigen(F);
1081 mini::SVector<ScalarType, Dims * Dims> gradPsiF;
1082 mini::SMatrix<ScalarType, Dims * Dims, Dims * Dims> hessPsiF;
1083 auto psiF = Psi.evalWithGradAndHessian(vecF, mug(g), lambdag(g), gradPsiF, hessPsiF);
1084 Eigen::Matrix<ScalarType, Dims * Dims, Dims * Dims> hessPsiFCorrected;
1085 math::linalg::FilterEigenvalues(
1086 ToEigen(hessPsiF),
1087 ToEigenvalueFilter(eSpdCorrection),
1088 hessPsiFCorrected);
1089 auto const GP = FromEigen(GPeg);
1090 auto GPsix = GradientWrtDofs<TElement, Dims>(gradPsiF, GP);
1091 auto HPsix = HessianWrtDofs<TElement, Dims>(FromEigen(hessPsiFCorrected), GP);
1092 auto heg = Hg.template block<kDofsPerElement, kDofsPerElement>(0, g * kDofsPerElement);
1093 Ug(g) += wg(g) * psiF;
1094 Gg.col(g) += wg(g) * ToEigen(GPsix);
1095 heg += wg(g) * ToEigen(HPsix);
1096 });
1097 }
1098 else if (eFlags == EElementElasticityComputationFlags::Gradient)
1099 {
1100 tbb::parallel_for(SizeType{0}, nQuadPts, [&](auto g) {
1101 auto const e = eg(g);
1102 auto const nodes = E.col(e);
1103 auto const xe = x.reshaped(Dims, nNodes)(Eigen::placeholders::all, nodes);
1104 auto const GPeg = GNeg.template block<kNodesPerElement, Dims>(0, g * Dims);
1105 Eigen::Matrix<ScalarType, Dims, Dims> const F = xe * GPeg;
1106 auto vecF = FromEigen(F);
1107 auto const gradPsiF = Psi.grad(vecF, mug(g), lambdag(g));
1108 auto const GP = FromEigen(GPeg);
1109 auto const GPsix = GradientWrtDofs<TElement, Dims>(gradPsiF, GP);
1110 Gg.col(g) += wg(g) * ToEigen(GPsix);
1111 });
1112 }
1113 else if (
1114 eFlags == (EElementElasticityComputationFlags::Gradient |
1115 EElementElasticityComputationFlags::Hessian))
1116 {
1117 tbb::parallel_for(SizeType{0}, nQuadPts, [&](auto g) {
1118 auto const e = eg(g);
1119 auto const nodes = E.col(e);
1120 auto const xe = x.reshaped(Dims, nNodes)(Eigen::placeholders::all, nodes);
1121 auto const GPeg = GNeg.template block<kNodesPerElement, Dims>(0, g * Dims);
1122 Eigen::Matrix<ScalarType, Dims, Dims> const F = xe * GPeg;
1123 auto vecF = FromEigen(F);
1124 mini::SVector<ScalarType, Dims * Dims> gradPsiF;
1125 mini::SMatrix<ScalarType, Dims * Dims, Dims * Dims> hessPsiF;
1126 Psi.evalWithGradAndHessian(vecF, mug(g), lambdag(g), gradPsiF, hessPsiF);
1127 Eigen::Matrix<ScalarType, Dims * Dims, Dims * Dims> hessPsiFCorrected;
1128 math::linalg::FilterEigenvalues(
1129 ToEigen(hessPsiF),
1130 ToEigenvalueFilter(eSpdCorrection),
1131 hessPsiFCorrected);
1132 auto const GP = FromEigen(GPeg);
1133 auto GPsix = GradientWrtDofs<TElement, Dims>(gradPsiF, GP);
1134 auto HPsix = HessianWrtDofs<TElement, Dims>(FromEigen(hessPsiFCorrected), GP);
1135 auto heg = Hg.template block<kDofsPerElement, kDofsPerElement>(0, g * kDofsPerElement);
1136 Gg.col(g) += wg(g) * ToEigen(GPsix);
1137 heg += wg(g) * ToEigen(HPsix);
1138 });
1139 }
1140 else /* if (eFlags == EElementElasticityComputationFlags::Hessian) */
1141 {
1142 tbb::parallel_for(SizeType{0}, nQuadPts, [&](auto g) {
1143 auto const e = eg(g);
1144 auto const nodes = E.col(e);
1145 Eigen::Matrix<ScalarType, Dims, kNodesPerElement> const xe =
1146 x.reshaped(Dims, nNodes)(Eigen::placeholders::all, nodes);
1147 Eigen::Matrix<ScalarType, kNodesPerElement, Dims> const GPeg =
1148 GNeg.template block<kNodesPerElement, Dims>(0, g * Dims);
1149 Eigen::Matrix<ScalarType, Dims, Dims> const F = xe * GPeg;
1150 auto vecF = FromEigen(F);
1151 auto hessPsiF = Psi.hessian(vecF, mug(g), lambdag(g));
1152 Eigen::Matrix<ScalarType, Dims * Dims, Dims * Dims> hessPsiFCorrected;
1153 math::linalg::FilterEigenvalues(
1154 ToEigen(hessPsiF),
1155 ToEigenvalueFilter(eSpdCorrection),
1156 hessPsiFCorrected);
1157 auto const GP = FromEigen(GPeg);
1158 mini::SMatrix<ScalarType, kDofsPerElement, kDofsPerElement> HPsix =
1159 HessianWrtDofs<TElement, Dims>(FromEigen(hessPsiFCorrected), GP);
1160 auto heg = Hg.template block<kDofsPerElement, kDofsPerElement>(0, g * kDofsPerElement);
1161 heg += wg(g) * ToEigen(HPsix);
1162 });
1163 }
1164}
1165
1166template <
1167 CElement TElement,
1168 int Dims,
1169 class TDerivedE,
1170 class TDerivedeg,
1171 class TDerivedHg,
1172 class TDerivedIn,
1173 class TDerivedOut>
1174inline void GemmHyperElastic(
1175 Eigen::DenseBase<TDerivedE> const& E,
1176 Eigen::Index nNodes,
1177 Eigen::DenseBase<TDerivedeg> const& eg,
1178 Eigen::MatrixBase<TDerivedHg> const& Hg,
1179 Eigen::MatrixBase<TDerivedIn> const& X,
1180 Eigen::DenseBase<TDerivedOut>& Y)
1181{
1182 PBAT_PROFILE_NAMED_SCOPE("pbat.fem.GemmHyperElastic");
1183 auto const numberOfDofs = Dims * nNodes;
1184 if (X.rows() != numberOfDofs or Y.rows() != numberOfDofs or X.cols() != Y.cols())
1185 {
1186 std::string const what = fmt::format(
1187 "Expected inputs and outputs to have rows |#nodes*kDims|={} and same number of "
1188 "columns, but got dimensions "
1189 "X,Y=({},{}), ({},{})",
1190 numberOfDofs,
1191 X.rows(),
1192 X.cols(),
1193 Y.rows(),
1194 Y.cols());
1195 throw std::invalid_argument(what);
1196 }
1197
1198 auto constexpr kDofsPerElement = Dims * TElement::kNodes;
1199 auto const nQuadPts = eg.size();
1200 for (auto c = 0; c < X.cols(); ++c)
1201 {
1202 for (auto g = 0; g < nQuadPts; ++g)
1203 {
1204 auto const e = eg(g);
1205 auto const nodes = E.col(e);
1206 auto const heg =
1207 Hg.template block<kDofsPerElement, kDofsPerElement>(0, g * kDofsPerElement);
1208 auto const xe = X.col(c).reshaped(Dims, nNodes)(Eigen::placeholders::all, nodes);
1209 auto ye = Y.col(c).reshaped(Dims, nNodes)(Eigen::placeholders::all, nodes);
1210 ye.reshaped() += heg * xe.reshaped();
1211 }
1212 }
1213}
1214
1215template <
1216 CElement TElement,
1217 int Dims,
1218 Eigen::StorageOptions Options,
1219 class TDerivedE,
1220 class TDerivedeg>
1221auto ElasticHessianSparsity(
1222 Eigen::DenseBase<TDerivedE> const& E,
1223 Eigen::Index nNodes,
1224 Eigen::DenseBase<TDerivedeg> const& eg)
1225 -> math::linalg::SparsityPattern<typename TDerivedE::Scalar, Options>
1226{
1227 PBAT_PROFILE_NAMED_SCOPE("pbat.fem.ElasticHessianSparsity");
1228 using IndexType = typename TDerivedE::Scalar;
1229 auto const numberOfQuadraturePoints = eg.size();
1230 auto const kNodesPerElement = TElement::kNodes;
1231 auto const kDofsPerElement = kNodesPerElement * Dims;
1232 auto const nDofs = Dims * nNodes;
1233 std::vector<IndexType> nonZeroRowIndices{};
1234 std::vector<IndexType> nonZeroColIndices{};
1235 nonZeroRowIndices.reserve(
1236 static_cast<std::size_t>(kDofsPerElement * kDofsPerElement * numberOfQuadraturePoints));
1237 nonZeroColIndices.reserve(
1238 static_cast<std::size_t>(kDofsPerElement * kDofsPerElement * numberOfQuadraturePoints));
1239 // Insert non-zero indices in the storage order of our Hg matrix of element hessians at
1240 // quadrature points
1241 for (auto g = 0; g < numberOfQuadraturePoints; ++g)
1242 {
1243 auto const e = eg(g);
1244 auto const nodes = E.col(e);
1245 for (auto j = 0; j < kNodesPerElement; ++j)
1246 {
1247 for (auto dj = 0; dj < Dims; ++dj)
1248 {
1249 for (auto i = 0; i < kNodesPerElement; ++i)
1250 {
1251 for (auto di = 0; di < Dims; ++di)
1252 {
1253 nonZeroRowIndices.push_back(Dims * nodes(i) + di);
1254 nonZeroColIndices.push_back(Dims * nodes(j) + dj);
1255 }
1256 }
1257 }
1258 }
1259 }
1260 math::linalg::SparsityPattern<IndexType, Options> GH;
1261 GH.Compute(nDofs, nDofs, nonZeroRowIndices, nonZeroColIndices);
1262 return GH;
1263}
1264
1265template <
1266 CElement TElement,
1267 int Dims,
1268 Eigen::StorageOptions Options,
1269 class TDerivedE,
1270 class TDerivedeg,
1271 class TDerivedHg>
1272auto HyperElasticHessian(
1273 Eigen::DenseBase<TDerivedE> const& E,
1274 Eigen::Index nNodes,
1275 Eigen::DenseBase<TDerivedeg> const& eg,
1276 Eigen::MatrixBase<TDerivedHg> const& Hg)
1277 -> Eigen::SparseMatrix<typename TDerivedHg::Scalar, Options, typename TDerivedE::Scalar>
1278{
1279 PBAT_PROFILE_NAMED_SCOPE("pbat.fem.HyperElasticHessian");
1280 using ScalarType = typename TDerivedHg::Scalar;
1281 using IndexType = typename TDerivedE::Scalar;
1282 using SparseMatrixType = Eigen::SparseMatrix<ScalarType, Options, IndexType>;
1283 using Triplet = Eigen::Triplet<ScalarType, IndexType>;
1284
1285 auto constexpr kNodesPerElement = TElement::kNodes;
1286 auto constexpr kDofsPerElement = kNodesPerElement * Dims;
1287 auto const nQuadPts = eg.size();
1288
1289 std::vector<Triplet> triplets{};
1290 triplets.reserve(static_cast<std::size_t>(kDofsPerElement * kDofsPerElement * nQuadPts));
1291
1292 auto const numberOfDofs = Dims * nNodes;
1293 SparseMatrixType H(numberOfDofs, numberOfDofs);
1294
1295 for (auto g = 0; g < nQuadPts; ++g)
1296 {
1297 auto const e = eg(g);
1298 auto const nodes = E.col(e);
1299 auto const heg =
1300 Hg.template block<kDofsPerElement, kDofsPerElement>(0, g * kDofsPerElement);
1301 for (auto j = 0; j < kNodesPerElement; ++j)
1302 {
1303 for (auto dj = 0; dj < Dims; ++dj)
1304 {
1305 for (auto i = 0; i < kNodesPerElement; ++i)
1306 {
1307 for (auto di = 0; di < Dims; ++di)
1308 {
1309 auto const ni = static_cast<IndexType>(Dims * nodes(i) + di);
1310 auto const nj = static_cast<IndexType>(Dims * nodes(j) + dj);
1311 triplets.emplace_back(ni, nj, heg(Dims * i + di, Dims * j + dj));
1312 }
1313 }
1314 }
1315 }
1316 }
1317 H.setFromTriplets(triplets.begin(), triplets.end());
1318 return H;
1319}
1320
1321template <
1322 CElement TElement,
1323 int Dims,
1324 class TDerivedE,
1325 class TDerivedHg,
1326 Eigen::StorageOptions Options,
1327 class TDerivedH>
1328void ToHyperElasticHessian(
1329 Eigen::DenseBase<TDerivedE> const& E,
1330 Eigen::Index nNodes,
1331 Eigen::DenseBase<TDerivedHg> const& Hg,
1332 math::linalg::SparsityPattern<typename TDerivedE::Scalar, Options> const& sparsity,
1333 Eigen::SparseCompressedBase<TDerivedH>& H)
1334{
1335 PBAT_PROFILE_NAMED_SCOPE("pbat.fem.ToHyperElasticHessian");
1336 using SpanType = std::span<Scalar const>;
1337 using SizeType = typename SpanType::size_type;
1338 sparsity.To(SpanType(Hg.data(), static_cast<SizeType>(Hg.size())), H);
1339}
1340
1341template <
1342 CElement TElement,
1343 int Dims,
1344 Eigen::StorageOptions Options,
1345 class TDerivedE,
1346 class TDerivedHg>
1347auto HyperElasticHessian(
1348 Eigen::DenseBase<TDerivedE> const& E,
1349 Eigen::Index nNodes,
1350 Eigen::DenseBase<TDerivedHg> const& Hg,
1351 math::linalg::SparsityPattern<typename TDerivedE::Scalar, Options> const& sparsity)
1352 -> Eigen::SparseMatrix<typename TDerivedHg::Scalar, Options, typename TDerivedE::Scalar>
1353{
1354 using SpanType = std::span<Scalar const>;
1355 using SizeType = typename SpanType::size_type;
1356 return sparsity.ToMatrix(SpanType(Hg.data(), static_cast<SizeType>(Hg.size())));
1357}
1358
1359template <
1360 CElement TElement,
1361 int Dims,
1362 class TDerivedE,
1363 class TDerivedeg,
1364 class TDerivedGg,
1365 class TDerivedOut>
1366inline void ToHyperElasticGradient(
1367 Eigen::DenseBase<TDerivedE> const& E,
1368 Eigen::Index nNodes,
1369 Eigen::DenseBase<TDerivedeg> const& eg,
1370 Eigen::MatrixBase<TDerivedGg> const& Gg,
1371 Eigen::PlainObjectBase<TDerivedOut>& G)
1372{
1373 PBAT_PROFILE_NAMED_SCOPE("pbat.fem.ToHyperElasticGradient");
1374 G.resize(Dims * nNodes, 1);
1375 G.setZero();
1376 auto constexpr kNodesPerElement = TElement::kNodes;
1377 auto const nQuadPts = eg.size();
1378 auto unvecG = G.reshaped(Dims, nNodes);
1379 for (auto g = 0; g < nQuadPts; ++g)
1380 {
1381 auto const e = eg(g);
1382 auto const nodes = E.col(e);
1383 auto const geg = Gg.col(g).reshaped(Dims, kNodesPerElement);
1384 auto gi = unvecG(Eigen::placeholders::all, nodes);
1385 gi += geg;
1386 }
1387}
1388
1389} // namespace pbat::fem
1390
1391#endif // PBAT_FEM_HYPERELASTICPOTENTIAL_H
DeformationGradient.h
Functions to compute deformation gradient and its derivatives.
pbat::math::linalg::SparsityPattern
Sparsity pattern precomputer to accelerate sparse matrix assembly.
Definition SparsityPattern.h:32
Eigen.h
Eigen adaptors for ranges.
pbat::fem::CElement
Reference finite element.
Definition Concepts.h:63
pbat::physics::CHyperElasticEnergy
Concept for hyperelastic energy.
Definition HyperElasticity.h:53
pbat::fem
Finite Element Method (FEM)
Definition Concepts.h:19
pbat::fem::ToHyperElasticHessian
void ToHyperElasticHessian(Eigen::DenseBase< TDerivedE > const &E, Eigen::Index nNodes, Eigen::DenseBase< TDerivedHg > const &Hg, math::linalg::SparsityPattern< typename TDerivedE::Scalar, Options > const &sparsity, Eigen::SparseCompressedBase< TDerivedH > &H)
Construct the hessian matrix using mesh with sparsity pattern.
Definition HyperElasticPotential.h:1328
pbat::fem::GemmHyperElastic
void GemmHyperElastic(Eigen::DenseBase< TDerivedE > const &E, Eigen::Index nNodes, Eigen::DenseBase< TDerivedeg > const &eg, Eigen::MatrixBase< TDerivedHg > const &Hg, Eigen::MatrixBase< TDerivedIn > const &X, Eigen::DenseBase< TDerivedOut > &Y)
Apply the hessian matrix as a linear operator .
Definition HyperElasticPotential.h:1174
pbat::fem::EElementElasticityComputationFlags
EElementElasticityComputationFlags
Bit-flag enum for element elasticity computation flags.
Definition HyperElasticPotential.h:54
pbat::fem::HessianWrtDofs
auto HessianWrtDofs(TMatrixHF const &HF, TMatrixGP const &GP) -> math::linalg::mini::SMatrix< ScalarType, TElement::kNodes *Dims, TElement::kNodes *Dims >
Computes hessian w.r.t. FEM degrees of freedom of scalar function , where is the jacobian of ,...
Definition DeformationGradient.h:231
pbat::fem::ToHyperElasticGradient
void ToHyperElasticGradient(Eigen::DenseBase< TDerivedE > const &E, Eigen::Index nNodes, Eigen::DenseBase< TDerivedeg > const &eg, Eigen::MatrixBase< TDerivedGg > const &Gg, Eigen::PlainObjectBase< TDerivedG > &G)
Compute the gradient vector into existing output.
pbat::fem::HyperElasticGradient
auto HyperElasticGradient(Eigen::DenseBase< TDerivedE > const &E, Eigen::Index nNodes, Eigen::DenseBase< TDerivedeg > const &eg, Eigen::MatrixBase< TDerivedGg > const &Gg) -> Eigen::Vector< typename TDerivedGg::Scalar, Eigen::Dynamic >
Compute the gradient vector (allocates output)
Definition HyperElasticPotential.h:651
pbat::fem::EHyperElasticSpdCorrection
EHyperElasticSpdCorrection
Bit-flag enum for SPD projection type.
Definition HyperElasticPotential.h:38
pbat::fem::GradientWrtDofs
auto GradientWrtDofs(TMatrixGF const &GF, TMatrixGP const &GP) -> math::linalg::mini::SVector< ScalarType, TElement::kNodes *Dims >
Computes gradient w.r.t. FEM degrees of freedom of scalar function , where is the jacobian of ,...
Definition DeformationGradient.h:126
pbat::fem::ToElementElasticity
void ToElementElasticity(Eigen::DenseBase< TDerivedE > const &E, Eigen::Index nNodes, Eigen::DenseBase< TDerivedeg > const &eg, Eigen::DenseBase< TDerivedwg > const &wg, Eigen::MatrixBase< TDerivedGNeg > const &GNeg, Eigen::DenseBase< TDerivedmug > const &mug, Eigen::DenseBase< TDerivedlambdag > const &lambdag, Eigen::MatrixBase< TDerivedx > const &x, Eigen::PlainObjectBase< TDerivedUg > &Ug, Eigen::PlainObjectBase< TDerivedGg > &Gg, Eigen::PlainObjectBase< TDerivedHg > &Hg, int eFlags=EElementElasticityComputationFlags::Potential, EHyperElasticSpdCorrection eSpdCorrection=EHyperElasticSpdCorrection::None)
Compute element elasticity and its derivatives at the given shape.
Definition HyperElasticPotential.h:934
pbat::fem::HyperElasticPotential
auto HyperElasticPotential(Eigen::DenseBase< TDerivedUg > const &Ug) -> typename TDerivedUg::Scalar
Compute the total elastic potential.
Definition HyperElasticPotential.h:852
pbat::fem::ToEigenvalueFilter
math::linalg::EEigenvalueFilter ToEigenvalueFilter(EHyperElasticSpdCorrection mode)
Convert EHyperElasticSpdCorrection to EEigenvalueFilter.
Definition HyperElasticPotential.cpp:5
pbat::fem::HyperElasticHessian
auto HyperElasticHessian(Eigen::DenseBase< TDerivedE > const &E, Eigen::Index nNodes, Eigen::DenseBase< TDerivedeg > const &eg, Eigen::MatrixBase< TDerivedHg > const &Hg) -> Eigen::SparseMatrix< typename TDerivedHg::Scalar, Options, typename TDerivedE::Scalar >
Construct the hessian matrix's sparse representation.
Definition HyperElasticPotential.h:1272
pbat::fem::ElasticHessianSparsity
auto ElasticHessianSparsity(Eigen::DenseBase< TDerivedE > const &E, Eigen::Index nNodes, Eigen::DenseBase< TDerivedeg > const &eg) -> math::linalg::SparsityPattern< typename TDerivedE::Scalar, Options >
Construct the hessian matrix's sparsity pattern.
Definition HyperElasticPotential.h:1221
pbat::math::linalg::mini
Mini linear algebra related functionality.
Definition Assign.h:12
pbat::math::linalg::EEigenvalueFilter
EEigenvalueFilter
Bit-flag enum for SPD projection type.
Definition FilterEigenvalues.h:12
pbat::math::linalg::FilterEigenvalues
bool FilterEigenvalues(Eigen::MatrixBase< TDerivedA > const &A, EEigenvalueFilter mode, Eigen::MatrixBase< TDerivedB > &B)
Filter eigenvalues of a symmetric matrix A and store the result in B.
Definition FilterEigenvalues.h:30
Profiling.h
Profiling utilities for the Physics-Based Animation Toolkit (PBAT)