Remove duplicated line and arc curves used in tests

Author: stfnp

solver/virtualbow_fem/src/elements/beam/geometry.rs

@@ -1,4 +1,4 @@
-use nalgebra::{SMatrix, SVector, matrix, vector};
+use nalgebra::{SMatrix, SVector};
 
 // This module contains common traits for describing a beam's geometry
 
@@ -57,124 +57,4 @@ pub trait CrossSection {
     // Returns the stress recovery matrices for the cross-section at relative position n and for implementation-specific points of interest.
     // When multiplied with the strain vector [epsilon, gamma, kappa], each matrix produces the normal stress at that point.
     fn stress_recovery(&self, n: f64) -> Vec<SVector<f64, 3>>;
-}
-
-// TODO: Duplicate with testutils/curves.rs?
-// Implementation of a linearly varying rectangular cross-section for use in tests
-pub struct RectangularSection {
-    pub w0: f64,
-    pub h0: f64,
-    pub w1: f64,
-    pub h1: f64,
-    pub ρ: f64,
-    pub E: f64,
-    pub G: f64,
-}
-
-impl CrossSection for RectangularSection {
-    fn stiffness(&self, n: f64) -> SMatrix<f64, 3, 3> {
-        let w = self.width(n);
-        let h = self.height(n);
-
-        let EA = self.E*w*h;
-        let GA = self.G*w*h;
-        let EI = self.E*w*h.powi(3)/12.0;
-
-        matrix![
-            EA, 0.0, 0.0;
-            0.0, GA, 0.0;
-            0.0, 0.0, EI;
-        ]
-    }
-
-    fn mass(&self, n: f64) -> SMatrix<f64, 3, 3> {
-        let w = self.width(n);
-        let h = self.height(n);
-
-        let ρA = self.ρ*w*h;
-        let ρI = self.ρ*w*h.powi(3)/12.0;
-
-        matrix![
-            ρA, 0.0, 0.0;
-            0.0, ρA, 0.0;
-            0.0, 0.0, ρI;
-        ]
-    }
-
-    fn width(&self, p: f64) -> f64 {
-        self.w0 + p*(self.w1 - self.w0)
-    }
-
-    fn height(&self, p: f64) -> f64 {
-        self.h0 + p*(self.h1 - self.h0)
-    }
-
-    fn strain_recovery(&self, _p: f64) -> Vec<SVector<f64, 3>> {
-        todo!();
-    }
-
-    fn stress_recovery(&self, _p: f64) -> Vec<SVector<f64, 3>> {
-        todo!();
-    }
-}
-
-// TODO: Duplicate with testutils/curves.rs?
-// Implementation of a straight line curve for use in tests
-pub struct LineCurve {
-    pub x: f64,
-    pub y: f64,
-    pub φ: f64,
-    pub l: f64,
-}
-
-impl PlanarCurve for LineCurve {
-    fn length(&self) -> f64 {
-        self.l
-    }
-
-    fn position(&self, s: f64) -> SVector<f64, 2> {
-        vector![
-            self.x + s*f64::cos(self.φ),
-            self.y + s*f64::sin(self.φ),
-        ]
-    }
-
-    fn angle(&self, _s: f64) -> f64 {
-        self.φ
-    }
-
-    fn curvature(&self, _s: f64) -> f64 {
-        0.0
-    }
-}
-
-// TODO: Duplicate with testutils/curves.rs?
-// Implementation of a circular arc curve for use in tests
-pub struct ArcCurve {
-    pub x: f64,
-    pub y: f64,
-    pub φ: f64,
-    pub l: f64,
-    pub r: f64,
-}
-
-impl PlanarCurve for ArcCurve {
-    fn length(&self) -> f64 {
-        self.l
-    }
-
-    fn position(&self, s: f64) -> SVector<f64, 2> {
-        vector![
-            self.x + self.r*(f64::sin(s/self.r + self.φ) - f64::sin(self.φ)),
-            self.y + self.r*(f64::cos(self.φ) - f64::cos(s/self.r + self.φ))
-        ]
-    }
-
-    fn angle(&self, s: f64) -> f64 {
-        self.φ + s/self.r
-    }
-
-    fn curvature(&self, _s: f64) -> f64 {
-        1.0/self.r
-    }
 }

solver/virtualbow_fem/src/elements/beam/linear.rs

@@ -196,26 +196,28 @@ mod tests {
     use approx::{assert_abs_diff_eq, assert_relative_eq};
     use iter_num_tools::lin_space;
     use nalgebra::{DMatrix, matrix, SMatrix, stack};
-    use crate::elements::beam::geometry::{ArcCurve, CrossSection, LineCurve, PlanarCurve, RectangularSection};
+    use crate::elements::beam::geometry::{CrossSection, PlanarCurve};
     use crate::elements::beam::linear::LinearBeamSegment;
+    use crate::testutils::curves::{ArcCurve, LineCurve};
+    use crate::testutils::sections::RectangularSection;
 
     #[test]
     fn test_linear_stiffness_matrix_straight() {
         // Computes the stiffness matrix of a straight beam segment with constant cross-section
         // and an arbitrary starting point and -angle and compares it to the exact solution
         // as well as the fem approximation.
 
-        let angle = 0.1;
         let length = 0.8;
-        let curve = LineCurve { x: 1.5, y: 2.0, φ: angle, l: length };
-        let section = RectangularSection { w0: 0.01, h0: 0.01, w1: 0.01, h1: 0.01, ρ: 7850.0, E: 210e9, G: 80e9 };
+        let angle = 0.1;
+        let curve = LineCurve::new([1.5, 2.0, angle], length);
+        let section = RectangularSection::new(7850.0, 210e9, 80e9, &[0.01], &[0.01], &[0.0]);
 
         let n_elements = 100;
         let segment_fem = LinearBeamSegmentFEM::new(&curve, &section, 0.0, curve.length(), n_elements);
         let segment_num = LinearBeamSegment::new(&curve, &section, 0.0, curve.length(), &segment_fem.s_eval);
 
         let C_sec = section.stiffness(0.0);
-        let K_ref = LinearBeamSegmentFEM::element_stiffness_matrix(C_sec[(0, 0)], C_sec[(1, 1)], C_sec[(2, 2)], curve.length(), angle);
+        let K_ref = LinearBeamSegmentFEM::element_stiffness_matrix(C_sec[(0, 0)], C_sec[(1, 1)], C_sec[(2, 2)], length, angle);
 
         // Check stiffness matrices
         assert_relative_eq!(segment_num.K, K_ref, max_relative=1e-9);
@@ -227,7 +229,7 @@ mod tests {
         }
 
         // Check total segment mass
-        assert_abs_diff_eq!(segment_num.m, section.ρ*section.w0*section.h0*length, epsilon=1e-12);
+        assert_abs_diff_eq!(segment_num.m, section.ρ*section.w[0]*section.h[0]*length, epsilon=1e-12);
     }
 
     #[test]
@@ -236,8 +238,8 @@ mod tests {
         // and an arbitrary starting point and -angle and compares it to the fem approximation.
 
         let length = 0.8;
-        let curve = ArcCurve { x: 1.5, y: 2.0, φ: 0.1, l: length, r: 0.4 };
-        let section = RectangularSection { w0: 0.01, h0: 0.01, w1: 0.01, h1: 0.01, ρ: 7850.0, E: 210e9, G: 80e9 };
+        let curve = ArcCurve::new([1.5, 2.0, 0.1], length, 0.4);
+        let section = RectangularSection::new(7850.0, 210e9, 80e9, &[0.01], &[0.01], &[0.0]);
 
         let n_elements = 100;
         let segment_fem = LinearBeamSegmentFEM::new(&curve, &section, 0.0, curve.length(), n_elements);
@@ -252,7 +254,7 @@ mod tests {
         }
 
         // Check total segment mass
-        assert_abs_diff_eq!(segment_num.m, section.ρ*section.w0*section.h0*length, epsilon=1e-12);
+        assert_abs_diff_eq!(segment_num.m, section.ρ*section.w[0]*section.h[0]*length, epsilon=1e-12);
     }
 
     #[test]
@@ -261,8 +263,8 @@ mod tests {
         // and an arbitrary starting point and -angle and compares it to the fem approximation.
 
         let length = 0.8;
-        let curve = ArcCurve { x: 1.5, y: 2.0, φ: 0.1, l: length, r: 0.4 };
-        let section = RectangularSection { w0: 0.01, h0: 0.01, w1: 0.005, h1: 0.005, ρ: 7850.0, E: 210e9, G: 80e9 };
+        let curve = ArcCurve::new([1.5, 2.0, 0.1], length, 0.4);
+        let section = RectangularSection::new(7850.0, 210e9, 80e9, &[0.01, 0.005], &[0.01, 0.005], &[0.0, 0.0]);
 
         let n_elements = 100;
         let segment_fem = LinearBeamSegmentFEM::new(&curve, &section, 0.0, curve.length(), n_elements);
@@ -278,8 +280,8 @@ mod tests {
 
         // Check total segment mass
         // Analytical volume from truncated pyramid: https://de.wikipedia.org/wiki/Pyramidenstumpf
-        let A0 = section.w0*section.h0;
-        let A1 = section.w1*section.h1;
+        let A0 = section.w[0]*section.h[0];
+        let A1 = section.w[1]*section.h[1];
         assert_abs_diff_eq!(segment_num.m, section.ρ*length/3.0*(A0 + f64::sqrt(A0*A1) + A1), epsilon=1e-12);
     }
 

solver/virtualbow_fem/src/testutils/curves.rs

@@ -3,66 +3,81 @@ use std::{f64::consts::FRAC_PI_2, f64::consts::PI};
 use nalgebra::{vector, SVector};
 use crate::elements::beam::geometry::PlanarCurve;
 
-// Simple straight line of given length for use in tests
-pub struct Line {
-    l: f64
+// Implementation of a straight line curve for use in tests
+pub struct LineCurve {
+    x: f64,
+    y: f64,
+    φ: f64,
+    l: f64,
 }
 
-impl Line {
-    pub fn new(l: f64) -> Self {
+impl LineCurve {
+    pub fn new(start: [f64; 3], l: f64) -> Self {
         Self {
-            l
+            x: start[0],
+            y: start[1],
+            φ: start[2],
+            l,
         }
     }
 }
 
-impl PlanarCurve for Line {
+impl PlanarCurve for LineCurve {
     fn length(&self) -> f64 {
         self.l
     }
 
     fn position(&self, s: f64) -> SVector<f64, 2> {
-        vector![s, 0.0]
+        vector![
+            self.x + s*f64::cos(self.φ),
+            self.y + s*f64::sin(self.φ),
+        ]
     }
 
     fn angle(&self, _s: f64) -> f64 {
-        0.0
+        self.φ
     }
 
     fn curvature(&self, _s: f64) -> f64 {
         0.0
     }
 }
 
-pub struct Arc {
+// Implementation of a circular arc curve for use in tests
+pub struct ArcCurve {
+    x: f64,
+    y: f64,
+    φ: f64,
     l: f64,
-    r: f64
+    r: f64,
 }
 
-impl Arc {
-    pub fn new(l: f64, r: f64) -> Self {
+impl ArcCurve {
+    pub fn new(start: [f64; 3], l: f64, r: f64) -> Self {
         Self {
+            x: start[0],
+            y: start[1],
+            φ: start[2],
             l,
             r
         }
     }
 }
 
-impl PlanarCurve for Arc {
+impl PlanarCurve for ArcCurve {
     fn length(&self) -> f64 {
         self.l
     }
 
     fn position(&self, s: f64) -> SVector<f64, 2> {
-        let φ = self.angle(s);
         vector![
-            self.r*f64::sin(φ),
-            self.r*(1.0 - f64::cos(φ))
+            self.x + self.r*(f64::sin(s/self.r + self.φ) - f64::sin(self.φ)),
+            self.y + self.r*(f64::cos(self.φ) - f64::cos(s/self.r + self.φ))
         ]
     }
 
     fn angle(&self, s: f64) -> f64 {
-        s/self.r
+        self.φ + s/self.r
     }
 
     fn curvature(&self, _s: f64) -> f64 {
@@ -89,8 +104,8 @@ impl Wave {
     fn center(&self, s: f64) -> (f64, f64, f64) {
         let n = s/self.l;
 
-        let i = f64::floor(n*(self.k as f64));            // Index of the current arc
-        let c = self.r*(1.0 + 2.0*i);                     // Center of the current arc
+        let i = f64::floor(n*(self.k as f64));        // Index of the current arc
+        let c = self.r*(1.0 + 2.0*i);                 // Center of the current arc
         let α = (s/self.l*(self.k as f64) - i)*PI;    // Angle from arc center
 
         (i, c, α)

solver/virtualbow_fem/src/testutils/sections.rs

@@ -3,16 +3,18 @@ use crate::elements::beam::geometry::CrossSection;
 
 // Simple rectangular cross-section for use in tests
 
-pub struct Section {
-    ρ: f64,
-    E: f64,
-    G: f64,
-    w: Vec<f64>,
-    h: Vec<f64>,
-    y: Vec<f64>
+pub struct RectangularSection {
+    pub ρ: f64,
+    pub E: f64,
+    pub G: f64,
+    pub w: Vec<f64>,
+    pub h: Vec<f64>,
+    pub y: Vec<f64>
 }
 
-impl Section {
+impl RectangularSection {
+    // TODO: &[(width, height, offset)]?
+    // TODO: More constructors (constant, linear, cubic)
     pub fn new(ρ: f64, E: f64, G: f64, w: &[f64], h: &[f64], y: &[f64]) -> Self {
         Self {
             ρ,
@@ -51,7 +53,7 @@ impl Section {
     }
 }
 
-impl CrossSection for Section {
+impl CrossSection for RectangularSection {
     fn stiffness(&self, n: f64) -> SMatrix<f64, 3, 3> {
         let w = self.width(n);
         let h = self.height(n);

solver/virtualbow_fem/tests/test_01_element_properties.rs

@@ -5,8 +5,9 @@ use virtualbow_fem::elements::beam::beam::BeamElement;
 use virtualbow_fem::elements::beam::linear::LinearBeamSegment;
 use virtualbow_fem::elements::string::StringElement;
 use virtualbow_fem::testutils::syschecks::assert_system_invariants;
-use virtualbow_fem::elements::beam::geometry::{ArcCurve, RectangularSection};
 use virtualbow_fem::system::dof::DofType;
+use virtualbow_fem::testutils::curves::ArcCurve;
+use virtualbow_fem::testutils::sections::RectangularSection;
 
 // These tests perform basic consistency checks on the various elements
 // See the utils::checks::check_system_invariants function for the details
@@ -88,8 +89,9 @@ fn beam_element() {
     let l = 0.6;
     let r = 0.4;
 
-    let curve = ArcCurve { x: 0.0, y: 0.0, φ: 0.0, l, r, };
-    let section = RectangularSection { w0: 0.01, h0: 0.01, w1: 0.005, h1: 0.005, ρ: 740.0, E: 11670e6, G: 8000e6 };
+    let curve = ArcCurve::new([0.0; 3], l, r);
+    let section = RectangularSection::new(740.0, 11670e6, 8000e6, &[0.01, 0.005], &[0.01, 0.005], &[0.0, 0.0]);
+
     let segment = LinearBeamSegment::new(&curve, &section, 0.0, l, &[]);
 
     let mut element = BeamElement::new(&segment);