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
ValanisLandelQuasiNewton.h
Go to the documentation of this file.
1
11
12#ifndef PBAT_SIM_ALGORITHM_PD_VALANISLANDELQUASINEWTON_H
13#define PBAT_SIM_ALGORITHM_PD_VALANISLANDELQUASINEWTON_H
14
15#include "pbat/common/Concepts.h"
16#include "pbat/physics/Enums.h"
17
18#include <Eigen/Core>
19#include <cmath>
20
21namespace pbat::sim::algorithm::pd {
22
36template <common::CArithmetic TScalar>
37TScalar ValanisLandelQuasiNewtonStiffness(
38 TScalar mu,
39 TScalar lambda,
40 Eigen::Vector<TScalar, 3> const& sigma,
41 TScalar sigmalo,
42 TScalar sigmahi,
43 physics::EHyperElasticEnergy energy)
44{
45 TScalar kstar{1};
46 switch (energy)
47 {
48 case physics::EHyperElasticEnergy::SaintVenantKirchhoff: {
49 TScalar const a0 = ((sigmahi) * (sigmahi) * (sigmahi));
50 TScalar const a1 = ((sigmahi) * (sigmahi));
51 TScalar const a2 = ((sigmalo) * (sigmalo));
52 TScalar const a3 = ((sigmalo) * (sigmalo) * (sigmalo));
53 TScalar const a4 = 12 * sigmahi;
54 TScalar const a5 = lambda * mu;
55 TScalar const a6 = a4 * a5;
56 TScalar const a7 = 12 * sigmalo;
57 TScalar const a8 = a5 * a7;
58 TScalar const a9 = ((lambda) * (lambda));
59 TScalar const a10 = ((mu) * (mu));
60 TScalar const a11 = 6 * a1;
61 TScalar const a12 = a11 * a5;
62 TScalar const a13 = 4 * a5;
63 TScalar const a14 = ((sigmahi) * (sigmahi) * (sigmahi) * (sigmahi));
64 TScalar const a15 = 3 * a5;
65 TScalar const a16 = 6 * a2;
66 TScalar const a17 = a16 * a5;
67 TScalar const a18 = ((sigmalo) * (sigmalo) * (sigmalo) * (sigmalo));
68 TScalar const a19 = ((sigma[1]) * (sigma[1]));
69 TScalar const a20 = ((sigma[2]) * (sigma[2]));
70 TScalar const a21 = sigma[1] * sigma[2];
71 TScalar const a22 = a21 * a5;
72 TScalar const a23 = 12 * a22;
73 TScalar const a24 = 8 * a22;
74 TScalar const a25 = a21 * a9;
75 TScalar const a26 = 12 * a25;
76 TScalar const a27 = 8 * a25;
77 TScalar const a28 = a19 * a20 * a9;
78 TScalar const a29 = 4 * a28;
79 TScalar const a30 = 3 * a28;
80 kstar = (1.0 / 8.0) *
81 std::abs(
82 (a0 * a13 + a0 * a24 + a0 * a27 + a0 * a29 - a1 * a23 - a1 * a26 -
83 a10 * a11 + a10 * a16 + a10 * a4 - a10 * a7 - a11 * a9 - a12 * a19 -
84 a12 * a20 + a12 - a13 * a3 - a14 * a15 - a14 * a30 + a15 * a18 + a16 * a9 +
85 a17 * a19 + a17 * a20 - a17 + a18 * a30 + a19 * a6 - a19 * a8 + a2 * a23 +
86 a2 * a26 + a20 * a6 - a20 * a8 - a24 * a3 - a27 * a3 - a29 * a3 + a4 * a9 -
87 a6 - a7 * a9 + a8) /
88 (lambda * (a0 - 3 * a1 + 3 * a2 - a3 + 3 * sigmahi - 3 * sigmalo)));
89 break;
90 }
91 case physics::EHyperElasticEnergy::StableNeoHookean: {
92 TScalar const a0 = ((sigmahi) * (sigmahi) * (sigmahi));
93 TScalar const a1 = ((sigmahi) * (sigmahi));
94 TScalar const a2 = ((sigmalo) * (sigmalo));
95 TScalar const a3 = ((sigmalo) * (sigmalo) * (sigmalo));
96 TScalar const a4 = 270 * lambda;
97 TScalar const a5 = 180 * mu;
98 TScalar const a6 = 135 * lambda;
99 TScalar const a7 = 60 * lambda;
100 TScalar const a8 = ((sigmahi) * (sigmahi) * (sigmahi) * (sigmahi));
101 TScalar const a9 = 45 * lambda;
102 TScalar const a10 = std::pow(sigmahi, 5);
103 TScalar const a11 = 6 * lambda;
104 TScalar const a12 = std::pow(sigmahi, 6);
105 TScalar const a13 = 5 * lambda;
106 TScalar const a14 = ((sigmalo) * (sigmalo) * (sigmalo) * (sigmalo));
107 TScalar const a15 = std::pow(sigmalo, 5);
108 TScalar const a16 = std::pow(sigmalo, 6);
109 TScalar const a17 = 90 * mu;
110 TScalar const a18 = 40 * mu;
111 TScalar const a19 = 30 * mu;
112 TScalar const a20 = 12 * mu;
113 TScalar const a21 = 10 * mu;
114 TScalar const a22 = ((sigma[1]) * (sigma[1]));
115 TScalar const a23 = 180 * lambda;
116 TScalar const a24 = a23 * sigmahi;
117 TScalar const a25 = ((sigma[1]) * (sigma[1]) * (sigma[1]) * (sigma[1]));
118 TScalar const a26 = a25 * sigmahi;
119 TScalar const a27 = 30 * lambda;
120 TScalar const a28 = ((sigma[2]) * (sigma[2]));
121 TScalar const a29 = ((sigma[2]) * (sigma[2]) * (sigma[2]) * (sigma[2]));
122 TScalar const a30 = a29 * sigmahi;
123 TScalar const a31 = a22 * sigmalo;
124 TScalar const a32 = a27 * sigmalo;
125 TScalar const a33 = a28 * sigmalo;
126 TScalar const a34 = 120 * mu;
127 TScalar const a35 = a34 * sigmahi;
128 TScalar const a36 = 60 * mu;
129 TScalar const a37 = a36 * sigmalo;
130 TScalar const a38 = a1 * a22;
131 TScalar const a39 = 90 * lambda;
132 TScalar const a40 = 15 * lambda;
133 TScalar const a41 = a1 * a40;
134 TScalar const a42 = a1 * a28;
135 TScalar const a43 = a22 * lambda;
136 TScalar const a44 = 20 * a0;
137 TScalar const a45 = a28 * lambda;
138 TScalar const a46 = 15 * a8;
139 TScalar const a47 = a2 * a39;
140 TScalar const a48 = a2 * a40;
141 TScalar const a49 = 20 * a3;
142 TScalar const a50 = 15 * a14;
143 TScalar const a51 = a1 * a19;
144 TScalar const a52 = a2 * a36;
145 TScalar const a53 = a19 * a2;
146 TScalar const a54 = a22 * a28;
147 kstar = (1.0 / 20.0) *
148 std::abs(
149 (a0 * a18 + a0 * a7 + a1 * a17 + a1 * a6 - a10 * a11 - a10 * a20 +
150 a11 * a15 + a12 * a13 + a12 * a21 - a13 * a16 + a14 * a19 + a14 * a9 +
151 a15 * a20 - a16 * a21 - a17 * a2 - a18 * a3 - a19 * a8 - a2 * a27 * a54 -
152 a2 * a6 + a22 * a24 + a22 * a35 + a22 * a47 + a22 * a52 - a23 * a31 -
153 a23 * a33 + a24 * a28 + a25 * a32 + a25 * a37 + a25 * a41 - a25 * a48 +
154 a25 * a51 - a25 * a53 - a26 * a27 - a26 * a36 + a27 * a28 * a38 -
155 a27 * a30 + a28 * a31 * a7 + a28 * a35 + a28 * a47 + a28 * a52 +
156 a29 * a32 + a29 * a37 + a29 * a41 - a29 * a48 + a29 * a51 - a29 * a53 -
157 a3 * a7 - a30 * a36 - a31 * a34 - a33 * a34 - a36 * a38 - a36 * a42 -
158 a38 * a39 - a39 * a42 - a4 * sigmahi + a4 * sigmalo - a43 * a44 +
159 a43 * a46 + a43 * a49 - a43 * a50 - a44 * a45 + a45 * a46 + a45 * a49 -
160 a45 * a50 - a5 * sigmahi + a5 * sigmalo - a54 * a7 * sigmahi - a8 * a9) /
161 (a0 - 3 * a1 + 3 * a2 - a3 + 3 * sigmahi - 3 * sigmalo));
162 break;
163 }
164 default: break;
165 }
166 return kstar;
167}
168
185template <
186 class TDerivedMu,
187 class TDerivedLambda,
188 class TDerivedSigma,
189 class TDerivedSigmaLo,
190 class TDerivedSigmaHi,
191 class TDerivedK>
192void ValanisLandelQuasiNewtonStiffness(
193 Eigen::DenseBase<TDerivedMu> const& mu,
194 Eigen::DenseBase<TDerivedLambda> const& lambda,
195 Eigen::DenseBase<TDerivedSigma> const& sigma,
196 Eigen::DenseBase<TDerivedSigmaLo> const& sigmalo,
197 Eigen::DenseBase<TDerivedSigmaHi> const& sigmahi,
198 physics::EHyperElasticEnergy energy,
199 Eigen::DenseBase<TDerivedK>& k)
200{
201 using ScalarType = typename TDerivedMu::Scalar;
202 k.resize(mu.rows(), mu.cols());
203 auto n = mu.size();
204 for (Eigen::Index i = 0; i < n; ++i)
205 {
206 k(i) = ValanisLandelQuasiNewtonStiffness<ScalarType>(
207 mu(i),
208 lambda(i),
209 sigma.col(i),
210 sigmalo(i),
211 sigmahi(i),
212 energy);
213 }
214}
215
216} // namespace pbat::sim::algorithm::pd
217
218#endif // PBAT_SIM_ALGORITHM_PD_VALANISLANDELQUASINEWTON_H
Concepts.h
Concepts for common types.
pbat::sim::algorithm::pd
Namespace for Projective Dynamics simulation algorithms.
pbat::sim::algorithm::pd::ValanisLandelQuasiNewtonStiffness
TScalar ValanisLandelQuasiNewtonStiffness(TScalar mu, TScalar lambda, Eigen::Vector< TScalar, 3 > const &sigma, TScalar sigmalo, TScalar sigmahi, physics::EHyperElasticEnergy energy)
Computes the element constraint stiffness for the quasi-Newton Projective Dynamics hessian approximat...
Definition ValanisLandelQuasiNewton.h:37