Skip to content

Commit 810f1eb

Browse files
jparr721jparr721
committed
working on straight strand assumption
1 parent ef53818 commit 810f1eb

10 files changed

+310
-85
lines changed

‎.idea/inspectionProfiles/Project_Default.xml

+1-1
Some generated files are not rendered by default. Learn more about customizing how changed files appear on GitHub.

‎CMakeLists.txt

+1-1
Original file line numberDiff line numberDiff line change
@@ -21,7 +21,7 @@ FetchContent_Declare(
2121
)
2222
FetchContent_MakeAvailable(libigl)
2323

24-
option(BUILD_TESTING "Build tests" OFF)
24+
option(BUILD_TESTING "Build tests" ON)
2525
if (BUILD_TESTING)
2626
enable_testing()
2727
include(GoogleTest)

‎Lib/Energy/Strand/DiscreteElasticRod.cc

+147-65
Original file line numberDiff line numberDiff line change
@@ -10,6 +10,7 @@ DiscreteElasticRod::DiscreteElasticRod(Mat<Real> vertices)
1010

1111
// Initialize the mass matrix - uniform mass right now.
1212
Vec<Real> masses = Vec<Real>::Ones(DOFs());
13+
masses.segment<3>(3) = Vec3<Real>::Constant(5);
1314
mInv = ConstructDiagonalSparseMatrix(masses);
1415

1516
Initialize();
@@ -19,16 +20,16 @@ void DiscreteElasticRod::Initialize() {
1920
// initial velocity of centerline
2021
velocities = Vec<Real>::Zero(vertices.rows());
2122

22-
// Thetas are per edge
23-
thetas = Vec<Real>::Zero(nRods);
23+
// Thetas are per bend frame
24+
thetas = Vec<Real>::Zero(nRods - 1);
2425

2526
// Lengths are per edge
2627
lengths = Vec<Real>::Zero(nRods);
2728
restLengths = Vec<Real>::Zero(nRods);
2829

29-
// Material frames are per edge
30-
m1s.resize(nRods);
31-
m2s.resize(nRods);
30+
// Material frames are per bend frame
31+
m1s.resize(nRods - 1);
32+
m2s.resize(nRods - 1);
3233

3334
for (int ii = 0; ii < nRods; ++ii) {
3435
const Real length = (vertices.row(ii + 1) - vertices.row(ii)).norm();
@@ -66,14 +67,16 @@ void DiscreteElasticRod::Initialize() {
6667
UpdateMaterialFrames();
6768

6869
// Init the bends
69-
bends.resize(nRods - 1);
70-
restBends.resize(nRods - 1);
70+
// bends.resize(nRods - 1);
71+
// restBends.resize(nRods - 1);
7172

72-
// Get the curvature update
73+
// Get the curvature update - don't bother with the quasistatic update here
74+
// since it's just 0 anyway.
7375
UpdateMaterialCurvatures();
76+
restCurvature = curvature;
7477

7578
// Set rest curvatures
76-
restBends = bends;
79+
// restBends = bends;
7780

7881
// Update the material curvature gradients
7982
UpdateKbGradients();
@@ -82,28 +85,31 @@ void DiscreteElasticRod::Initialize() {
8285
UpdateHolonomyGradient();
8386
}
8487

85-
auto DiscreteElasticRod::ComputeCenterlineForces() -> Vec<Real> {
88+
auto DiscreteElasticRod::ComputeCenterlineForcesGeneral() -> Vec<Real> {
8689
Vec<Real> R = Vec<Real>::Zero(DOFs());
8790

8891
// J is a 90 degree rotation matrix. Always
8992
Mat2<Real> J;
9093
J << 0, -1, 1, 0;
9194

9295
// Compute the derivative of the curvature
93-
const Mat2<Real> Bhat = Mat2<Real>::Identity() * BENDING_MODULUS;
96+
const Mat2<Real> Bhat = Mat2<Real>::Identity() * gBendingModulus;
9497

9598
// Formula in 7.1, the General case - dEdx
9699
for (int ii = 0; ii < nRods - 1; ++ii) {
97100
Vec3<Real> vertexForce = Vec3<Real>::Zero();
98101

99102
for (int jj = ii; jj <= ii + 1; ++jj) {
100-
const auto &w =
101-
jj == ii ? bends.at(ii).prevCurvature : bends.at(ii).nextCurvature;
102-
const auto &wbar = jj == ii ? restBends.at(ii).prevCurvature
103-
: bends.at(ii).nextCurvature;
103+
const auto &w = curvature.at(ii);
104+
const auto &wbar = restCurvature.at(ii);
105+
// const auto &w =
106+
// jj == ii ? bends.at(ii).prevCurvature :
107+
// bends.at(ii).nextCurvature;
108+
// const auto &wbar = jj == ii ? restBends.at(ii).prevCurvature
109+
// : bends.at(ii).nextCurvature;
104110

105-
const Mat2x3<Real> gradW = ComputeGradW(gradientkbs.at(ii), m1s.at(jj),
106-
m2s.at(jj), holonomy.at(ii), w);
111+
const Mat2x3<Real> gradW = ComputeGradOmega(
112+
gradientkbs.at(ii), m1s.at(jj), m2s.at(jj), holonomy.at(ii), w);
107113

108114
vertexForce += gradW.transpose() * Bhat * (w - wbar);
109115
}
@@ -116,6 +122,16 @@ auto DiscreteElasticRod::ComputeCenterlineForces() -> Vec<Real> {
116122
return R;
117123
}
118124

125+
auto DiscreteElasticRod::ComputeCenterlineForcesStraight() -> Vec<Real> {
126+
Vec<Real> R = Vec<Real>::Zero(DOFs());
127+
128+
for (int ii = 0; ii < nRods - 2; ++ii) {
129+
R.segment<3>(3 * ii + 1) += ComputepEpxi(ii);
130+
}
131+
132+
return R;
133+
}
134+
119135
void DiscreteElasticRod::UpdateBishopFrames() {
120136
ASSERT2(!kbs.empty());
121137
for (int ii = 1; ii < nRods; ++ii) {
@@ -140,7 +156,7 @@ void DiscreteElasticRod::UpdateBishopFrames() {
140156
}
141157

142158
void DiscreteElasticRod::UpdateMaterialFrames() {
143-
for (int ii = 0; ii < nRods; ++ii) {
159+
for (int ii = 0; ii < nRods - 1; ++ii) {
144160
const Real theta = thetas(ii);
145161
const Real cosTheta = std::cos(theta);
146162
const Real sinTheta = std::sin(theta);
@@ -154,61 +170,50 @@ void DiscreteElasticRod::UpdateMaterialFrames() {
154170
}
155171

156172
void DiscreteElasticRod::UpdateMaterialCurvatures() {
173+
curvature.clear();
174+
curvature.resize(nRods - 1);
157175
for (int ii = 0; ii < nRods - 1; ++ii) {
158176
const Vec3<Real> kb = kbs.at(ii);
159177
const Vec3<Real> m1 = m1s.at(ii);
160178
const Vec3<Real> m2 = m2s.at(ii);
161179

162-
bends.at(ii).prevCurvature = ComputeW(kb, m1, m2);
163-
164-
if (ii > 0) {
165-
bends.at(ii - 1).nextCurvature = ComputeW(kb, m1, m2);
166-
}
180+
curvature.at(ii) = ComputeOmega(kb, m1, m2);
181+
//
182+
// bends.at(ii).prevCurvature = ComputeOmega(kb, m1, m2);
183+
//
184+
// if (ii > 0) {
185+
// bends.at(ii - 1).nextCurvature = ComputeOmega(kb, m1, m2);
186+
// }
167187
}
168188
}
169189

170190
void DiscreteElasticRod::UpdateQuasistaticMaterialFrame() {
171191
// J is a 90 degree rotation matrix
172192
Mat2<Real> J;
173193
J << 0, -1, 1, 0;
174-
const Mat2<Real> Bhat = Mat2<Real>::Identity() * BENDING_MODULUS;
194+
const Mat2<Real> Bhat = Mat2<Real>::Identity() * gBendingModulus;
175195

176-
Real dEdTheta = 0.0;
177-
Vec<Real> lowerDiagonal = Vec<Real>::Zero(nRods);
178-
Vec<Real> upperDiagonal = Vec<Real>::Zero(nRods);
179-
Vec<Real> diagonal = Vec<Real>::Zero(nRods);
196+
Vec<Real> dEdTheta = Vec<Real>::Zero(nRods - 1);
197+
Vec<Real> lowerDiagonal = Vec<Real>::Zero(nRods - 1);
198+
Vec<Real> upperDiagonal = Vec<Real>::Zero(nRods - 1);
199+
Vec<Real> diagonal = Vec<Real>::Zero(nRods - 1);
180200

181-
const auto computeGradWi = [&J, &Bhat](const Vec2<Real> &w,
182-
const Vec2<Real> &wbar,
183-
Real restLength) -> Real {
184-
const Real invRestLength = 1.0 / restLength;
185-
return invRestLength * w.transpose() * J * Bhat * (w - wbar);
186-
};
201+
// Run a single newton step
202+
const Vec<Real> twistForce = ComputeTwistingForce();
203+
ComputeTwistingHessian(upperDiagonal, diagonal, lowerDiagonal);
187204

188-
// 8 newton steps
189-
for (int iter = 0; iter < 8; ++iter) {
190-
for (int ii = 0; ii < nRods - 1; ++ii) {
191-
// Previous curvature W_j + 2beta * (mj/lj)
192-
{
193-
const Vec2<Real> w = bends.at(ii).prevCurvature;
194-
const Vec2<Real> wbar = restBends.at(ii).prevCurvature;
195-
const Real restLength = restLengths(ii);
196-
const Real mj = thetas(ii + 1) - thetas(ii);
197-
const Real Wi = computeGradWi(w, wbar, restLength);
198-
dEdTheta += Wi + 2 * BENDING_MODULUS * mj / restLength;
199-
}
200-
201-
// Next curvature W_j+1 = 2 * beta * (mj+1/lj+1)
202-
if (ii > 0) {
203-
const Vec2<Real> w = bends.at(ii - 1).nextCurvature;
204-
const Vec2<Real> wbar = restBends.at(ii - 1).nextCurvature;
205-
const Real restLength = restLengths(ii - 1);
206-
const Real mj = thetas(ii) - thetas(ii - 1);
207-
const Real Wi = computeGradWi(w, wbar, restLength);
208-
dEdTheta += Wi + 2 * BENDING_MODULUS * mj / restLength;
209-
}
210-
}
211-
}
205+
// Compute the theta update
206+
thetas =
207+
FactorTriDiagonalMatrix(upperDiagonal, diagonal, lowerDiagonal, dEdTheta);
208+
209+
std::cout << "new thetas: " << thetas.transpose() << std::endl;
210+
211+
#ifndef NDEBUG
212+
std::cout << "thetas.norm(): " << thetas.norm() << std::endl;
213+
#endif
214+
215+
// Compute the quasistatic frames with the new theta values
216+
UpdateMaterialFrames();
212217
}
213218

214219
void DiscreteElasticRod::UpdateKbGradients() {
@@ -264,16 +269,17 @@ void DiscreteElasticRod::Computekbs() {
264269
}
265270
}
266271

267-
auto DiscreteElasticRod::ComputeW(const Vec3<Real> &kb, const Vec3<Real> &m1,
268-
const Vec3<Real> &m2) -> Vec2<Real> {
272+
auto DiscreteElasticRod::ComputeOmega(const Vec3<Real> &kb,
273+
const Vec3<Real> &m1,
274+
const Vec3<Real> &m2) -> Vec2<Real> {
269275
return {kb.dot(m2), -kb.dot(m1)};
270276
}
271277

272-
auto DiscreteElasticRod::ComputeGradW(const Mat3<Real> &gradkb,
273-
const Vec3<Real> &m1,
274-
const Vec3<Real> &m2,
275-
const Vec3<Real> &gradpsi,
276-
const Vec2<Real> &w) -> Mat2x3<Real> {
278+
auto DiscreteElasticRod::ComputeGradOmega(const Mat3<Real> &gradkb,
279+
const Vec3<Real> &m1,
280+
const Vec3<Real> &m2,
281+
const Vec3<Real> &gradpsi,
282+
const Vec2<Real> &w) -> Mat2x3<Real> {
277283
// J is a 90 degree rotation
278284
Mat2<Real> J;
279285
J << 0, -1, 1, 0;
@@ -284,3 +290,79 @@ auto DiscreteElasticRod::ComputeGradW(const Mat3<Real> &gradkb,
284290

285291
return M * gradkb - J * w * gradpsi.transpose();
286292
}
293+
294+
auto DiscreteElasticRod::ComputeTwistingForce() -> Vec<Real> {
295+
Mat2<Real> J;
296+
J << 0, -1, 1, 0;
297+
const Mat2<Real> Bhat = Mat2<Real>::Identity() * gBendingModulus;
298+
299+
Vec<Real> forces = Vec<Real>::Zero(nRods - 1);
300+
301+
const auto gradW = [&J, &Bhat, this](int j) -> Real {
302+
const auto &w = curvature.at(j);
303+
const auto &wbar = restCurvature.at(j);
304+
const Real restLength = restLengths(j);
305+
const Real invRestLength = 1.0 / restLength;
306+
return invRestLength * w.transpose() * J * Bhat * (w - wbar);
307+
};
308+
309+
// Directly from equation 7.
310+
for (int j = 0; j < nRods - 2; ++j) {
311+
const auto &Wj = gradW(j);
312+
const auto &Wj1 = gradW(j + 1);
313+
const Real thetaDiff = thetas(j) - thetas(j + 1);
314+
315+
forces(j) = (Wj + Wj1) * 2 * gBendingModulus * thetaDiff / restLengths(j);
316+
}
317+
318+
#ifndef NDEBUG
319+
std::cout << "forces: " << forces.transpose() << std::endl;
320+
#endif
321+
322+
return forces;
323+
}
324+
325+
void DiscreteElasticRod::ComputeTwistingHessian(Vec<Real> &upper,
326+
Vec<Real> &center,
327+
Vec<Real> &lower) {
328+
Mat<Real> hessian;
329+
330+
Mat2<Real> J;
331+
J << 0, -1, 1, 0;
332+
const Mat2<Real> Bhat = Mat2<Real>::Identity() * gBendingModulus;
333+
const auto hessW = [&J, &Bhat, this](int j) -> Real {
334+
const auto &w = curvature.at(j);
335+
const auto &wbar = restCurvature.at(j);
336+
const Real restLength = restLengths(j);
337+
const Real invRestLength = 1.0 / restLength;
338+
339+
const Real lhs = invRestLength * w.transpose() * J.transpose() * Bhat * w;
340+
const Real rhs = invRestLength * w.transpose() * Bhat * (w - wbar);
341+
return lhs - rhs;
342+
};
343+
344+
// Components of the hessian tri-diagonal matrix
345+
lower = Vec<Real>::Zero(nRods - 1);
346+
upper = Vec<Real>::Zero(nRods - 1);
347+
center = Vec<Real>::Zero(nRods - 1);
348+
349+
for (int j = 0; j < nRods - 2; ++j) {
350+
lower(j) = -2 * gBendingModulus / restLengths(j);
351+
upper(j) = -2 * gBendingModulus / restLengths(j + 1);
352+
353+
center(j) =
354+
hessW(j) + hessW(j + 1) +
355+
2 * gBendingModulus * (1 / restLengths(j) + 1 / restLengths(j + 1));
356+
}
357+
}
358+
359+
auto DiscreteElasticRod::ComputepEpxi(int j) -> Vec3<Real> {
360+
const auto scaledTwistContrib = 2 * gTwistingModulus / restLengths(j);
361+
const auto &gradKb = gradientkbs.at(j);
362+
const auto &kb = kbs.at(j);
363+
const Real totalRestLength = restLengths.sum();
364+
const auto &gradHolo = holonomy.at(j);
365+
const Real thetaDiff = thetas(thetas.rows() - 1) - thetas(0);
366+
return scaledTwistContrib * gradKb.transpose() * kb +
367+
(gBendingModulus * thetaDiff) / totalRestLength * gradHolo;
368+
}

‎Lib/Energy/Strand/DiscreteElasticRod.h

+23-15
Original file line numberDiff line numberDiff line change
@@ -2,8 +2,8 @@
22

33
#include <LibMath.h>
44

5-
#define BENDING_MODULUS 1.0
6-
#define TWISTING_MODULUS 0.075;
5+
constexpr Real gBendingModulus = 1.0;
6+
constexpr Real gTwistingModulus = 0.075;
77

88
INLINE auto RotationMatrixAroundNormal(const Vec3<Real> &axisAngle)
99
-> Mat3<Real> {
@@ -82,14 +82,11 @@ class DiscreteElasticRod {
8282
std::vector<Mat3<Real>> gradientkbs;
8383

8484
// Keeps track of material curvature
85-
std::vector<Bend> bends;
86-
std::vector<Bend> restBends;
85+
// std::vector<Bend> bends;
86+
// std::vector<Bend> restBends;
8787

88-
// std::vector<Vec2<Real>> nextCurvature;
89-
// std::vector<Vec2<Real>> prevCurvature;
90-
//
91-
// std::vector<Vec2<Real>> restNextCurvature;
92-
// std::vector<Vec2<Real>> restPrevCurvature;
88+
std::vector<Vec2<Real>> curvature;
89+
std::vector<Vec2<Real>> restCurvature;
9390

9491
Vec<Real> velocities;
9592
Vec<Real> thetas;
@@ -109,7 +106,8 @@ class DiscreteElasticRod {
109106
}
110107

111108
void Initialize();
112-
auto ComputeCenterlineForces() -> Vec<Real>;
109+
auto ComputeCenterlineForcesGeneral() -> Vec<Real>;
110+
auto ComputeCenterlineForcesStraight() -> Vec<Real>;
113111

114112
/**
115113
* Update the bishop frame with a rotation matrix
@@ -123,9 +121,19 @@ class DiscreteElasticRod {
123121

124122
void Computekbs();
125123

126-
auto ComputeW(const Vec3<Real> &kb, const Vec3<Real> &m1,
127-
const Vec3<Real> &m2) -> Vec2<Real>;
128-
auto ComputeGradW(const Mat3<Real> &gradkb, const Vec3<Real> &m1,
129-
const Vec3<Real> &m2, const Vec3<Real> &gradpsi,
130-
const Vec2<Real> &w) -> Mat2x3<Real>;
124+
auto ComputeOmega(const Vec3<Real> &kb, const Vec3<Real> &m1,
125+
const Vec3<Real> &m2) -> Vec2<Real>;
126+
auto ComputeGradOmega(const Mat3<Real> &gradkb, const Vec3<Real> &m1,
127+
const Vec3<Real> &m2, const Vec3<Real> &gradpsi,
128+
const Vec2<Real> &w) -> Mat2x3<Real>;
129+
130+
/**
131+
* Partial of E wrt theta^j from equation 7.
132+
* @return
133+
*/
134+
auto ComputeTwistingForce() -> Vec<Real>;
135+
void ComputeTwistingHessian(Vec<Real> &upper, Vec<Real> &center,
136+
Vec<Real> &lower);
137+
138+
auto ComputepEpxi(int j) -> Vec3<Real>;
131139
};

0 commit comments

Comments
 (0)