Expose geometry computation via the solver api

Author: stfnp

gui/source/solver/API.cpp

@@ -53,6 +53,12 @@ void save_result(const BowResult& result, const std::string& path) {
     check_response(response);
 }
 
+DiscreteLimbGeometry compute_geometry(const BowModel& model) {
+    std::vector<uint8_t> data = json::to_msgpack(model);
+    Response response = ffi::compute_geometry(data.data(), data.size());
+    return parse_response<DiscreteLimbGeometry>(response);
+}
+
 BowResult simulate_model(const BowModel& model, Mode mode, bool (*callback)(Mode, double)) {
     std::vector<uint8_t> data = json::to_msgpack(model);
     Response response = ffi::simulate_model(data.data(), data.size(), mode, callback);

gui/source/solver/API.hpp

@@ -31,4 +31,6 @@ BowResult load_result(const std::string& path);
 
 void save_result(const BowResult& result, const std::string& path);
 
+DiscreteLimbGeometry compute_geometry(const BowModel& model);
+
 BowResult simulate_model(const BowModel& model, Mode mode, bool (*callback)(Mode, double));

gui/source/solver/BowResult.hpp

@@ -24,6 +24,19 @@ struct nlohmann::adl_serializer<std::optional<T>> {
     }
 };
 
+struct DiscreteLimbGeometry {
+    std::vector<double> n_eval;
+    std::vector<double> w_eval;
+    std::vector<double> h_eval;
+};
+
+NLOHMANN_DEFINE_TYPE_NON_INTRUSIVE(
+    DiscreteLimbGeometry,
+    n_eval,
+    w_eval,
+    h_eval
+)
+
 struct LayerInfo {
     std::string name;
 };

gui/source/tests/main.cpp

@@ -36,6 +36,11 @@ TEST_CASE("save-result-file") {
     REQUIRE_THROWS(save_result(result, TEST_DATA_DIR + "/temp/nonexistent/result.res"));
 }
 
+TEST_CASE("compute-geometry") {
+    BowModel model = new_model();
+    DiscreteLimbGeometry geometry = compute_geometry(model);
+}
+
 TEST_CASE("simulate-model") {
     BowModel model = new_model();
     BowResult result = simulate_model(model, Mode::Dynamic, [](Mode mode, double progress) {

rust/virtualbow/src/geometry.rs

@@ -1,6 +1,7 @@
 use iter_num_tools::lin_space;
 use itertools::Itertools;
 use nalgebra::{DMatrix, DVector, SVector, vector};
+use serde::{Deserialize, Serialize};
 use crate::errors::ModelError;
 use crate::input::{BowModel, HandleReference};
 use crate::profile::profile::{CurvePoint, ProfileCurve};
@@ -13,6 +14,26 @@ pub struct LimbGeometry {
     pub section: LayeredCrossSection,    // Limb cross sections
 }
 
+// TODO: Return values s_nodes, u_node might not be needed if the evaluation works properly
+#[derive(Serialize, Deserialize, Default, PartialEq, Debug, Clone)]
+pub struct DiscreteLimbGeometry {
+    pub segments: Vec<LinearBeamSegment>,    // Linear beam segment properties
+    pub n_nodes: Vec<f64>,                   // Relative lengths of the element nodes
+    pub s_nodes: Vec<f64>,                   // Arc lengths of the element nodes
+    pub p_nodes: Vec<SVector<f64, 3>>,       // Positions (x, y, φ) of the element nodes
+    pub y_nodes: Vec<DVector<f64>>,          // Layer bounds at nodes (y in cross section coordinates)
+
+    pub n_eval: Vec<f64>,                    // Relative lengths at which the limb quantities are evaluated
+    pub s_eval: Vec<f64>,                    // Arc lengths at which the limb quantities are evaluated
+    pub p_eval: Vec<SVector<f64, 3>>,        // Positions (x, y, φ) of the evaluation points
+    pub y_eval: Vec<DVector<f64>>,           // Layer bounds at eval points (y in cross section coordinates)
+    pub w_eval: Vec<f64>,                    // Widths at eval points
+    pub h_eval: Vec<f64>,                    // Total heights at eval points
+
+    pub strain_eval: Vec<DMatrix<f64>>,      // Strain evaluation matrices for each evaluation point
+    pub stress_eval: Vec<DMatrix<f64>>,      // Stress evaluation matrices for each evaluation point
+}
+
 impl LimbGeometry {
     pub fn new(input: &BowModel) -> Result<Self, ModelError> {
         // Section properties according to layers, materials and alignment to the profile curve
@@ -59,14 +80,14 @@ impl LimbGeometry {
     pub fn discretize(&self, n_eval_points: usize, n_elements: usize) -> DiscreteLimbGeometry {
         // Arc lengths and normalized positions along the profile where the element nodes are placed
         let s_nodes = lin_space(self.profile.s_start()..=self.profile.s_end(), n_elements + 1).collect_vec();
-        let p_nodes = s_nodes.iter().map(|&s| self.profile.normalize(s)).collect_vec();
-        let u_nodes = s_nodes.iter().map(|&s| self.profile.point(s)).collect_vec();
-        let y_nodes = p_nodes.iter().map(|&p| self.section.layer_bounds(p).0).collect_vec();
+        let n_nodes = s_nodes.iter().map(|&s| self.profile.normalize(s)).collect_vec();
+        let p_nodes = s_nodes.iter().map(|&s| self.profile.point(s)).collect_vec();
+        let y_nodes = n_nodes.iter().map(|&n| self.section.layer_bounds(n).0).collect_vec();
 
         // Equidistant evaluation points along the length of the limb
         let s_eval = lin_space(self.profile.s_start()..=self.profile.s_end(), n_eval_points).collect_vec();
-        let p_eval = s_eval.iter().map(|&s| self.profile.normalize(s)).collect_vec();
-        let y_eval = p_eval.iter().map(|&p| self.section.layer_bounds(p).0).collect_vec();
+        let n_eval = s_eval.iter().map(|&s| self.profile.normalize(s)).collect_vec();
+        let y_eval = n_eval.iter().map(|&n| self.section.layer_bounds(n).0).collect_vec();
 
         let segments = s_nodes.iter().tuple_windows().enumerate().map(|(i, (&s0, &s1))| {
             // TODO: Better solution for numerical issues?
@@ -86,45 +107,47 @@ impl LimbGeometry {
             LinearBeamSegment::new(&self.profile, &self.section, s0, s1, &s_eval)
         }).collect();
 
-        let strain_eval = p_eval.iter().map(|&p| self.section.strain_eval(p)).collect();
-        let stress_eval = p_eval.iter().map(|&p| self.section.stress_eval(p)).collect();
+        let p_eval = s_eval.iter().map(|&s| self.profile.point(s)).collect();
+        let w_eval = n_eval.iter().map(|&n| self.section.width(n)).collect();
+        let h_eval = n_eval.iter().map(|&n| self.section.height(n)).collect();
 
-        let position = s_eval.iter().map(|&s| self.profile.point(s)).collect();
-        let width = p_eval.iter().map(|&p| self.section.width(p)).collect();
-        let height = p_eval.iter().map(|&p| self.section.height(p)).collect();
+        let strain_eval = n_eval.iter().map(|&n| self.section.strain_eval(n)).collect();
+        let stress_eval = n_eval.iter().map(|&n| self.section.stress_eval(n)).collect();
 
         DiscreteLimbGeometry {
             segments,
+            n_nodes,
             s_nodes,
-            u_nodes,
+            p_nodes,
             y_nodes,
+            n_eval,
             s_eval,
             y_eval,
             strain_eval,
             stress_eval,
-            position,
-            width,
-            height
+            p_eval,
+            w_eval,
+            h_eval
         }
     }
 }
 
-// TODO: Return values s_nodes, u_node might not be needed if the evaluation works properly
-pub struct DiscreteLimbGeometry {
-    pub segments: Vec<LinearBeamSegment>,    // Linear beam segment properties
-    pub s_nodes: Vec<f64>,                   // Arc lengths of the element nodes
-    pub u_nodes: Vec<SVector<f64, 3>>,       // Positions (x, y, φ) of the element nodes
-    pub y_nodes: Vec<DVector<f64>>,          // Layer bounds at nodes (y in cross section coordinates)
+impl TryInto<Vec<u8>> for DiscreteLimbGeometry {
+    type Error = ModelError;
 
-    pub s_eval: Vec<f64>,                    // Arc lengths at which the limb quantities are evaluated (positions, forces, ...)
-    pub y_eval: Vec<DVector<f64>>,           // Layer bounds at eval points (y in cross section coordinates)
-    pub strain_eval: Vec<DMatrix<f64>>,      // Strain evaluation matrices for each evaluation point
-    pub stress_eval: Vec<DMatrix<f64>>,      // Stress evaluation matrices for each evaluation point
+    // Conversion into MsgPack byte array
+    fn try_into(self) -> Result<Vec<u8>, Self::Error> {
+        rmp_serde::to_vec_named(&self).map_err(ModelError::OutputEncodeMsgPackError)  // TODO: Bett error type?
+    }
+}
+
+impl TryFrom<&[u8]> for DiscreteLimbGeometry {
+    type Error = ModelError;
 
-    // TODO: Unify with rest
-    pub position: Vec<SVector<f64, 3>>,
-    pub width: Vec<f64>,
-    pub height: Vec<f64>
+    // Conversion from MsgPack byte array
+    fn try_from(value: &[u8]) -> Result<Self, Self::Error> {
+        rmp_serde::from_slice(value).map_err(ModelError::OutputDecodeMsgPackError)
+    }
 }
 
 #[cfg(test)]

rust/virtualbow/src/sections/section.rs

@@ -195,39 +195,40 @@ impl LayeredCrossSection {
         &self.layers
     }
 
-    // Computes the layer boundaries at arc length s. Also returns the heights as a byproduct.
-    pub fn layer_bounds(&self, p: f64) -> (DVector<f64>, DVector<f64>) {
-        let h = self.layer_heights(p);
+    // Computes the layer boundaries at normalized length n in [0, 1].
+    // Also returns the layer heights as a byproduct.
+    pub fn layer_bounds(&self, n: f64) -> (DVector<f64>, DVector<f64>) {
+        let h = self.layer_heights(n);
         (&self.stacking*&h, h)
     }
 
-    // Computes the section boundaries (belly, back) at arc length s.
+    // Computes the section boundaries (belly, back) at normalized length n.
     // Equivalent to first and last layer bound.
-    pub fn section_bounds(&self, p: f64) -> (f64, f64) {
-        let (y, _) = self.layer_bounds(p);
+    pub fn section_bounds(&self, n: f64) -> (f64, f64) {
+        let (y, _) = self.layer_bounds(n);
         (y[0], y[y.len() - 1])
     }
 
-    // Evaluates the heights of the individual layers at arc length s and returns them as a vector
-    pub fn layer_heights(&self, p: f64) -> DVector<f64> {
+    // Evaluates the heights of the individual layers at normalized length n and returns them as a vector
+    pub fn layer_heights(&self, n: f64) -> DVector<f64> {
         DVector::<f64>::from_fn(self.layers.len(), |i, _| {
-            self.layers[i].height.value(p, Extrapolation::Constant)
+            self.layers[i].height.value(n, Extrapolation::Constant)
         })
     }
 }
 
 impl CrossSection for LayeredCrossSection {
-    fn ρA(&self, p: f64) -> f64 {
-        let w = self.width.value(p, Extrapolation::Constant);
+    fn ρA(&self, n: f64) -> f64 {
+        let w = self.width.value(n, Extrapolation::Constant);
         self.layers.iter().map(|layer| {
-            let h = layer.height.value(p, Extrapolation::Constant);
+            let h = layer.height.value(n, Extrapolation::Constant);
             layer.material.density*w*h
         }).sum()
     }
 
-    fn ρI(&self, p: f64) -> f64 {
-        let w = self.width(p);
-        let (y, h) = self.layer_bounds(p);
+    fn ρI(&self, n: f64) -> f64 {
+        let w = self.width(n);
+        let (y, h) = self.layer_bounds(n);
 
         self.layers.iter().enumerate().map(|(i, layer)| {
             let A = w*h[i];
@@ -237,12 +238,12 @@ impl CrossSection for LayeredCrossSection {
         }).sum()
     }
 
-    fn C(&self, p: f64) -> SMatrix<f64, 3, 3> {
-        let w = self.width(p);
-        let (y, h) = self.layer_bounds(p);
+    fn C(&self, n: f64) -> SMatrix<f64, 3, 3> {
+        let w = self.width(n);
+        let (y, h) = self.layer_bounds(n);
 
         let Cee = self.layers.iter().map(|layer| {
-            let h = layer.height.value(p, Extrapolation::Constant);
+            let h = layer.height.value(n, Extrapolation::Constant);
             layer.material.youngs_modulus*w*h
         }).sum();
 
@@ -260,7 +261,7 @@ impl CrossSection for LayeredCrossSection {
         }).sum();
 
         let Cγγ = self.layers.iter().map(|layer| {
-            let h = layer.height.value(p, Extrapolation::Constant);
+            let h = layer.height.value(n, Extrapolation::Constant);
             layer.material.shear_modulus*w*h
         }).sum();
 
@@ -271,19 +272,19 @@ impl CrossSection for LayeredCrossSection {
         ]
     }
 
-    fn width(&self, p: f64) -> f64 {
-        self.width.value(p, Extrapolation::Constant)
+    fn width(&self, n: f64) -> f64 {
+        self.width.value(n, Extrapolation::Constant)
     }
 
     // Total height as the sum of all layers
-    fn height(&self, p: f64) -> f64 {
-        self.layers.iter().map(|layer| layer.height.value(p, Extrapolation::Constant)).sum()
+    fn height(&self, n: f64) -> f64 {
+        self.layers.iter().map(|layer| layer.height.value(n, Extrapolation::Constant)).sum()
     }
 
     // Strain evaluation matrix. Produces normal strains at back and belly of each layer.
-    fn strain_eval(&self, p: f64) -> DMatrix<f64> {
+    fn strain_eval(&self, n: f64) -> DMatrix<f64> {
         // Layer bounds are the points of interest
-        let (y, _) = self.layer_bounds(p);
+        let (y, _) = self.layer_bounds(n);
 
         // Normal strain is epsilon - kappa*y
         DMatrix::from_fn(2*self.layers.len(), 3, |i, j| {
@@ -297,9 +298,9 @@ impl CrossSection for LayeredCrossSection {
     }
 
     // Stress evaluation matrix. Produces stresses at back and belly of each layer.
-    fn stress_eval(&self, p: f64) -> DMatrix<f64> {
+    fn stress_eval(&self, n: f64) -> DMatrix<f64> {
         // Layer bounds are the points of interest
-        let (y, _) = self.layer_bounds(p);
+        let (y, _) = self.layer_bounds(n);
 
         // Normal stress is E*(epsilon - kappa*y)
         DMatrix::from_fn(2*self.layers.len(), 3, |i, j| {

rust/virtualbow/src/simulation.rs

@@ -77,7 +77,7 @@ impl<'a> Simulation<'a> {
 
         let mut system = System::new();
 
-        let limb_nodes: Vec<Node> = geometry.u_nodes.iter().enumerate().map(|(i, u)| {
+        let limb_nodes: Vec<Node> = geometry.p_nodes.iter().enumerate().map(|(i, u)| {
             system.create_node(u, &[i != 0; 3])    // First node is fixed, all others
         }).collect();
 
@@ -236,9 +236,9 @@ impl<'a> Simulation<'a> {
         let common = Common {
             limb: LimbInfo {
                 length: simulation.geometry.s_eval.clone(),
-                position: simulation.geometry.position.clone(),
-                width: simulation.geometry.width.clone(),
-                height: simulation.geometry.height.clone(),
+                position: simulation.geometry.p_eval.clone(),
+                width: simulation.geometry.w_eval.clone(),
+                height: simulation.geometry.h_eval.clone(),
                 bounds: simulation.geometry.y_eval.iter().map(|y| y.data.clone().into()).collect(),  // TODO: Uglyyy
             },
             layers,

rust/virtualbow_ffi/src/api.rs

@@ -1,5 +1,6 @@
 use std::path::Path;
 use virtualbow::errors::ModelError;
+use virtualbow::geometry::LimbGeometry;
 use virtualbow::input::BowModel;
 use virtualbow::input::BowModelVersion;
 use virtualbow::output::BowResult;
@@ -47,6 +48,15 @@ pub fn save_result<P>(data: &[u8], path: P) -> Result<(), String>
     result.save(path).map_err(|e| e.to_string())
 }
 
+pub fn compute_geometry(data: &[u8]) -> Result<Vec<u8>, String> {
+    let model = BowModel::try_from(data).map_err(|e| e.to_string())?;
+    let geometry = LimbGeometry::new(&model).map_err(|e| e.to_string())?;
+    let discretized = geometry.discretize(model.settings.num_limb_eval_points, model.settings.num_limb_elements);
+    let data = discretized.try_into().map_err(|e: ModelError| e.to_string())?;
+
+    Ok(data)
+}
+
 pub fn simulate_model<F>(data: &[u8], mode: SimulationMode, callback: F) -> Result<Vec<u8>, String>
     where F: Fn(SimulationMode, f64) -> bool
 {

rust/virtualbow_ffi/src/lib.rs

@@ -120,6 +120,17 @@ pub unsafe extern "C" fn save_result(data: *const u8, size: usize, path: *const
     }
 }
 
+#[no_mangle]
+pub unsafe extern "C" fn compute_geometry(data: *const u8, size: usize) -> Response {
+    let data = std::slice::from_raw_parts(data, size);
+    let result = api::compute_geometry(&data);
+
+    match result {
+        Ok(res) => Response::data(res),
+        Err(msg) => Response::error(msg),
+    }
+}
+
 #[repr(C)]
 pub enum Mode {
     Static,

rust/virtualbow_ffi/virtualbow.hpp

@@ -33,6 +33,8 @@ Response load_result(const char *path);
 
 Response save_result(const uint8_t *data, uintptr_t size, const char *path);
 
+Response compute_geometry(const uint8_t *data, uintptr_t size);
+
 Response simulate_model(const uint8_t *data, uintptr_t size, Mode mode, bool (*callback)(Mode, double));
 
 void free_response(Response response);