working strand bending!!!

Author: jparr721

.idea/misc.xml

@@ -38,7 +38,7 @@ macro(package_add_scene FILENAME)
 endmacro()
 package_add_scene(BunnyDropExplicit)
 package_add_scene(StrandDropExplicit)
-package_add_scene(StrandTwistExplicit)
+package_add_scene(StrandBendExplicit)
 
 add_executable(viewer main.cc)
 target_link_libraries(viewer PRIVATE ${OPENGL_LIBRARIES} GLUT::GLUT ${PROJECT_NAME})

CMakeLists.txt

@@ -2,7 +2,7 @@
 
 DiscreteElasticRod::DiscreteElasticRod(Mat<Real> vertices)
     : vertices(std::move(vertices)), nRods(vertices.rows() - 1),
-      massSpring(new MassSpring(100000, 0.25)) {
+      massSpring(new MassSpring(100000, 1)) {
   // At least 3 nodes is required
   ASSERT(this->vertices.rows() >= 3, "At least 3 nodes are required");
 
@@ -109,10 +109,14 @@ auto DiscreteElasticRod::ComputeCenterlineForcesGeneral() -> Vec<Real> {
 auto DiscreteElasticRod::ComputeCenterlineForcesStraight() -> Vec<Real> {
   Vec<Real> R = Vec<Real>::Zero(DOFs());
 
-//  for (int ii = 0; ii < nRods - 2; ++ii) {
-//    const Vec3<Real> vertexForce = ComputepEpxi(ii);
-//    R.segment<3>(3 * ii + 1) += vertexForce;
-//  }
+  for (int ii = 0; ii < nRods - 2; ++ii) {
+    const Vec3<Real> vertexForce = ComputepEpxi(ii);
+    R.segment<3>(3 * ii + 1) += vertexForce;
+  }
+
+#ifndef NDEBUG
+  std::cout << "forces: " << R.transpose() << std::endl;
+#endif
 
   // Compute spring-mass forces
   for (int ii = 0; ii < nRods; ++ii) {
@@ -286,10 +290,6 @@ auto DiscreteElasticRod::ComputeTwistingForce() -> Vec<Real> {
                  (thetaDiff1 / restEdges.at(j + 1).length));
   }
 
-#ifndef NDEBUG
-  std::cout << "forces: " << forces.transpose() << std::endl;
-#endif
-
   return forces;
 }
 
@@ -329,13 +329,18 @@ void DiscreteElasticRod::ComputeTwistingHessian(Vec<Real> &upper,
 }
 
 auto DiscreteElasticRod::ComputepEpxi(int j) -> Vec3<Real> {
-  const auto scaledTwistContrib = 2 * gTwistingModulus / restEdges.at(j).length;
+  Mat2<Real> J;
+  J << 0, -1, 1, 0;
+  const Mat2<Real> Bhat = Mat2<Real>::Identity() * gBendingModulus;
+
+  const auto &restLength = restEdges.at(j).length;
   const auto &gradKb = kbGradients.at(j);
-  const auto &kb = kbs.at(j);
-  const auto &gradHolo = holonomyGradients.at(j);
-  const Real thetaDiff = thetas(thetas.rows() - 1) - thetas(0);
-  return scaledTwistContrib * gradKb.transpose() * kb +
-         (gBendingModulus * thetaDiff) / totalRestLength * gradHolo;
+
+  const auto gradW = ComputeGradOmega(gradKb, m1s.at(j), m2s.at(j),
+                                      holonomyGradients.at(j), curvature.at(j));
+  const Vec2<Real> curvatureDiff = (curvature.at(j) - restCurvature.at(j));
+
+  return gradW.transpose() * Bhat * curvatureDiff / restLength;
 }
 
 auto DiscreteElasticRod::ComputepEpTheta() -> Real {

Lib/Energy/Strand/DiscreteElasticRod.cc

@@ -4,7 +4,7 @@
 #include <Energy/Spring/MassSpring.h>
 #include <LibMath.h>
 
-constexpr Real gBendingModulus = 0.1;
+constexpr Real gBendingModulus = 0.000001;
 constexpr Real gTwistingModulus = 0.7;
 
 class DiscreteElasticRod {

Lib/Energy/Strand/DiscreteElasticRod.h

@@ -84,7 +84,8 @@ class ForwardEulerStrand : public TimestepperStrand {
     this->strandMesh->der->UpdateLengths();
     this->strandMesh->der->Computekbs();
     this->strandMesh->der->UpdateBishopFrames();
-    this->strandMesh->der->UpdateQuasistaticMaterialFrame();
+    this->strandMesh->der->UpdateBishopFrames();
+    this->strandMesh->der->UpdateMaterialCurvatures();
     this->strandMesh->der->UpdateKbGradients();
     this->strandMesh->der->UpdateHolonomyGradient();
   }

Lib/ForwardEulerStrand.h

@@ -18,6 +18,11 @@ class StrandMesh {
   [[nodiscard]] INLINE auto DOFs() const -> int { return der->DOFs(); }
   INLINE void SetPositions(const Vec<Real> &x) { der->SetPositions(x); }
   INLINE auto Positions() const -> Vec<Real> { return der->Positions(); }
+  INLINE void TranslatePinnedX(Real dx, int ii) {
+    if (pinned[ii] == 1) {
+      der->vertices(ii, 0) = der->restVertices(ii, 0) + dx;
+    }
+  }
   void Draw();
 
   /**

Lib/StrandMesh.h

@@ -16,8 +16,8 @@ bool gCameraZooming = false;
 bool gCameraPanning = false;
 bool gAnimating = false;
 bool gSingleStep = false;
-bool gShowGrid = true;
-bool gShowCenterAxes = true;
+bool gShowGrid = false;
+bool gShowCenterAxes = false;
 bool gSaveFrame = false;
 
 int gSteps = 0;
@@ -182,26 +182,67 @@ void Display() {
   glFlush();
 }
 
+static Vec<Real> gTranslationStart = Vec<Real>::LinSpaced(5000, 0, 2.5);
+static Vec<Real> gTranslationEnd = Vec<Real>::LinSpaced(5000, -2.5, 0);
+static int gTranslationIndex = 0;
+static int gFramesSaved = 0;
+static int gStopFrame = 10000;
 static void GlutIdle() {
   if (gStopped) {
     return;
   }
 
   if (gAnimating || gSingleStep) {
-    Vec3<Real> gravity;
-    gravity.setZero();
+    Vec3<Real> gravity(0, -0.981, 0);
     ++gSteps;
 
-    gMesh->der->thetas(0) += 0.01;
-    gMesh->der->thetas(gMesh->der->thetas.rows() - 1) += 0.01;
+    gIntegrator->AddGravity(gravity);
     gIntegrator->Step();
 
     if (gSingleStep) {
       gSingleStep = !gSingleStep;
     }
+
+    if (gSteps > 1000 && gTranslationIndex < gTranslationStart.rows()) {
+      gMesh->TranslatePinnedX(gTranslationStart(gTranslationIndex), 0);
+      gMesh->TranslatePinnedX(
+          gTranslationEnd.reverse().eval()(gTranslationIndex), 4);
+      ++gTranslationIndex;
+    }
+  }
+
+  if (gSteps == gStopFrame) {
+    gAnimating = false;
+    gStopped = true;
   }
 
   if (gSteps % 10 == 0 && gAnimating) {
+    if (gSaveFrame) {
+
+
+      char filename[512];
+      sprintf(filename, "/Users/jarredparr/Downloads/output/frame_%05i.obj",
+              gFramesSaved);
+      ++gFramesSaved;
+
+      std::ofstream file(filename);
+
+      const auto &vertices = gMesh->der->vertices;
+      for (int jj = 0; jj < vertices.rows(); ++jj) {
+        // Write vertices.row(ii) to the file
+        file << "v " << vertices.row(jj) << std::endl;
+      }
+
+      file << "l ";
+      for (int jj = 0; jj < vertices.rows(); ++jj) {
+        int index = (jj + 1);
+        file << index << " ";
+      }
+
+      file << std::endl;
+      file.close();
+    }
+
     glutPostRedisplay();
   }
 }
@@ -222,12 +263,14 @@ auto main(int argc, char **argv) -> int {
   // Set the mesh
   Mat<Real> v(5, 3);
   for (int ii = 0; ii < v.rows(); ++ii) {
-    v.row(ii) = Vec3<Real>(ii, 2, 0);
+    v.row(ii) = Vec3<Real>(ii, 0, 0);
   }
 
   gMesh = std::make_shared<StrandMesh>(v);
-  gIntegrator =
-      std::make_shared<ForwardEulerStrand>(gMesh, nullptr, 1 / 100'000, 5, 0.3);
+  gMesh->pinned(0) = 1;
+  gMesh->pinned(v.rows() - 1) = 1;
+  gIntegrator = std::make_shared<ForwardEulerStrand>(gMesh, nullptr,
+                                                     1.0 / 1'000.0, 5, 0.3);
 
   glutInit(&argc, argv);
   glutInitDisplayMode(GLUT_DEPTH | GLUT_RGBA);

Scenes/StrandTwistExplicit.cc Scenes/StrandBendExplicit.cc

@@ -6,14 +6,14 @@
 
 app = typer.Typer()
 
-scene = """
+scene_drop = """
 <scene version="3.0.0">
 
     <!-- Defaults, these can be set via the command line: -Darg=value -->
 
     <default name="spp" value="300"/>
-    <default name="resx" value="1920"/>
-    <default name="resy" value="1080"/>
+    <default name="resx" value="1000"/>
+    <default name="resy" value="800"/>
     <default name="max_depth" value="8"/>
 
 
@@ -156,11 +156,11 @@ def make_render_scene(folder: str, radius: float):
 
     for file in track(os.listdir(folder)):
         # Convert frames to splines in the designated folder
-        if file.startswith("frame"):
+        if file.startswith("frame") and file.endswith(".obj"):
             load_obj_into_spline(os.path.join(folder, file), radius)
 
     # Now create the scene, first by copying the template
-    s = scene
+    s = scene_drop
 
     # Now create the spline entry and dump the scenes
     for file in track(os.listdir(folder)):
@@ -170,8 +170,9 @@ def make_render_scene(folder: str, radius: float):
             s += make_obj(os.path.join(folder, file))
 
         s += "</scene>\n"
-        print(s)
-        return
+        fn = os.path.splitext(file)[0]
+        with open(os.path.join(folder, fn + ".xml"), "w") as f:
+            f.write(s)
 
 
 if __name__ == "__main__":