core: add tension as optional fiber input

flatironinstitute · Feb 2, 2023 · 57d0a0d · 57d0a0d
1 parent 2fbf8fe
commit 57d0a0d
Show file tree

Hide file tree

Showing 4 changed files with 10 additions and 0 deletions.
diff --git a/include/fiber.hpp b/include/fiber.hpp
@@ -118,6 +118,9 @@ class Fiber {
             x_ = Eigen::MatrixXd::Zero(3, n_nodes_);
             x_.row(0) = Eigen::ArrayXd::LinSpaced(n_nodes_, 0, 1.0).transpose();
         }
+        if (tension_.size() != n_nodes_) {
+            tension_ = Eigen::VectorXd::Zero(n_nodes_);
+        }
 
         xs_.resize(3, n_nodes_);
         xss_.resize(3, n_nodes_);

diff --git a/src/core/fiber.cpp b/src/core/fiber.cpp
@@ -31,7 +31,10 @@ const std::string Fiber::BC_name[] = {"Force", "Torque", "Velocity", "AngularVel
 Fiber::Fiber(toml::value &fiber_table, double eta) {
     std::vector<double> x_array = toml::find<std::vector<double>>(fiber_table, "x");
     n_nodes_ = x_array.size() / 3;
+    std::vector<double> tension_array = toml::find_or(fiber_table, "tension", std::vector<double>{});
     x_ = Eigen::Map<Eigen::MatrixXd>(x_array.data(), 3, n_nodes_);
+    if (tension_array.size())
+        tension_ = Eigen::Map<Eigen::VectorXd>(tension_array.data(), n_nodes_);
 
     init();
 

diff --git a/src/core/system.cpp b/src/core/system.cpp
@@ -262,6 +262,7 @@ void prep_state_for_solver() {
         x_fibers.segment(offset + n * 0, n) = fiber.x_.row(0);
         x_fibers.segment(offset + n * 1, n) = fiber.x_.row(1);
         x_fibers.segment(offset + n * 2, n) = fiber.x_.row(2);
+        x_fibers.segment(offset + n * 3, n) = fiber.tension_;
         offset += 4 * n;
     }
     MatrixXd force_fibers = fc_.apply_fiber_force(x_fibers);

diff --git a/src/skelly_sim/skelly_config.py b/src/skelly_sim/skelly_config.py
@@ -273,6 +273,8 @@ class Fiber():
     x : List[float], default: :obj:`[]`, units: :obj:`μm`
         List of node positions in [x0,y0,z0,x1,y1,z1...] order. Extreme care must be taken when setting this since the length constraint
         can generate massive tensions with poor input. See examples.
+    tension : List[float], default: :obj:`[]`, units: :obj:`pN·μm⁻¹`
+        Optional list of node tensions [T_0,T_1,T_2,...T_{n_nodes-1}]
     """
     n_nodes: int = 32
     parent_site: int = -1
@@ -282,6 +284,7 @@ class Fiber():
     length: float = 1.0
     minus_clamped: bool = False
     x: List[float] = field(default_factory=list)
+    tension: List[float] = field(default_factory=list)
 
     def fill_node_positions(self, x0: np.array, normal: np.array):
         """