calculate 2d beam element stiffness matrix

src/bin/two_d.rs

@@ -0,0 +1,70 @@
+use fea::two_d::beam_elements::*;
+use fea::two_d::*;
+
+fn main() {
+    let student_id_last_digit = 9 as f32;
+    let d = student_id_last_digit / 100.0;
+    let r = d/2.0;
+    let e = 210e9;
+    let p = (25.0 - student_id_last_digit)*1000.0;
+    let theta = 30.0_f32.to_radians();
+    let l = 2.0/theta.cos();
+
+    let p1 = Point{
+        id: 1,
+        pos: Vector(0.0, 0.0),
+        bc: BoundaryCondition::Fixed,
+        beams: vec![],
+    };
+
+    let p2 = Point {
+        id: 2,
+        pos: Vector(0.0, l * theta.sin()),
+        bc: BoundaryCondition::Force(Vector(0.0, 1000.0)),
+        beams: vec![],
+    };
+
+    let l_e3 = 2.0/theta.sin();
+    println!("l_e3 {l_e3}");
+
+    let p3 = Point {
+        id: 3,
+        pos: Vector(0.0, (l * theta.sin()) + (l_e3 * theta.cos())),
+        bc: BoundaryCondition::Fixed,
+        beams: vec![],
+    };
+
+    let p4 = Point {
+        id: 4,
+        pos: Vector(l * theta.cos(), l * theta.sin()),
+        bc: BoundaryCondition::Force(Vector::from_mag_angle(p, 45_f32.to_radians())),
+        beams: vec![]
+    };
+
+    let mut w = World::from(vec![p1, p2, p3, p4]);
+
+    w.link(1,4, NewBeam {
+        cross_section: CrossSection::Circular(r),
+        material: fea::Material {
+            yield_stress: 95e6,
+            youngs_modulus: e,
+        }
+    }).unwrap();
+    w.link(2,4, NewBeam {
+        cross_section: CrossSection::Circular(r),
+        material: fea::Material {
+            yield_stress: 95e6,
+            youngs_modulus: e,
+        }
+    }).unwrap();
+    w.link(3,4, NewBeam {
+        cross_section: CrossSection::Circular(r),
+        material: fea::Material {
+            yield_stress: 95e6,
+            youngs_modulus: e,
+        }
+    }).unwrap();
+
+    println!("{:?}", w);
+    println!("{:?}", w.stiffness_matrix().unwrap());
+}

src/lib.rs

@@ -0,0 +1,7 @@
+pub mod two_d;
+
+#[derive(Debug)]
+pub struct Material {
+    pub youngs_modulus: f32,
+    pub yield_stress: f32,
+}

src/two_d.rs

@@ -0,0 +1,14 @@
+pub mod beam_elements;
+
+#[derive(Debug)]
+pub struct Vector(pub f32, pub f32);
+
+impl Vector {
+    pub fn distance(&self, other: &Self) -> f32 {
+        ((self.0 - other.0).powi(2) + (self.1 - other.1).powi(2)).sqrt()
+    }
+
+    pub fn from_mag_angle(magnitude: f32, angle: f32) -> Self {
+        Self (magnitude * angle.cos(), magnitude * angle.sin())
+    }
+}

src/two_d/beam_elements.rs

@@ -0,0 +1,25 @@
+mod beam;
+mod cross_section;
+mod world;
+
+pub use beam::*;
+pub use cross_section::*;
+pub use world::*;
+
+use super::Vector;
+use crate::Material;
+
+#[derive(Debug)]
+pub struct Point {
+    pub id: usize,
+    pub pos: Vector,
+    pub bc: BoundaryCondition,
+    pub beams: Vec<usize>,
+}
+
+#[derive(Debug)]
+pub enum BoundaryCondition {
+    Free,
+    Fixed,
+    Force(Vector),
+}

src/two_d/beam_elements/beam.rs

@@ -0,0 +1,68 @@
+use super::*;
+
+#[derive(Debug)]
+pub struct Beam{
+    pub id: usize,
+    pub points: (usize, usize),
+    pub material: Material,
+    pub cross_section: CrossSection,
+}
+
+pub struct NewBeam {
+    pub material: Material,
+    pub cross_section: CrossSection,
+}
+
+impl Beam {
+    pub fn new(p1: usize, p2: usize, props: NewBeam) -> Beam {
+        let b = Beam {
+            id: 0,
+            points: (p1, p2),
+            material: props.material,
+            cross_section: props.cross_section,
+        };
+
+        return b;
+    }
+
+    fn get_points<'a>(&self, world: &'a World) -> (&'a Point, &'a Point) {
+        let p1 = world.points.get(&self.points.0).unwrap();
+        let p2 = world.points.get(&self.points.1).unwrap();
+
+        (p1, p2)
+    }
+
+    pub fn stiffness(&self, world: &World) -> f32 {
+        let s = self.cross_section.area() * self.material.youngs_modulus / self.length(world);
+
+        println!("{:?} {}", self, s);
+
+        s
+
+    }
+
+    pub fn length(&self, world: &World) -> f32 {
+        let (p1, p2) = self.get_points(world);
+        let l = p1.pos.distance(&p2.pos);
+
+        println!("{:?} {}", self, l);
+
+        l
+    }
+
+    pub fn angle(&self, world: &World) -> f32 {
+        let (p1, p2) = self.get_points(world);
+        let dx = p1.pos.0 - p2.pos.0;
+        let dy = p1.pos.1 - p2.pos.1;
+
+        (dy/dx).atan()
+    }
+
+    pub fn other_point(&self, p: usize) -> usize {
+        if p == self.points.0{
+            return self.points.1;
+        } else {
+            return self.points.0;
+        }
+    }
+}

src/two_d/beam_elements/cross_section.rs

@@ -0,0 +1,21 @@
+use std::f32::consts::PI;
+
+#[derive(Debug)]
+pub enum CrossSection {
+    Area(f32),
+    Circular(f32),
+    Rectangular(f32, f32),
+}
+
+impl CrossSection {
+    pub fn area(&self) -> f32 {
+        let a = match self {
+            Self::Area(a) => *a,
+            Self::Circular(r) => PI * r.powi(2),
+            Self::Rectangular(l1, l2) => l1 * l2,
+        };
+        println!("{:?} {}", self, a);
+
+        a
+    }
+}

src/two_d/beam_elements/world.rs

@@ -0,0 +1,97 @@
+use super::*;
+use ndarray::prelude::*;
+use std::collections::HashMap;
+
+#[derive(Debug)]
+pub struct World {
+    pub points: HashMap<usize, Point>,
+    pub beams: HashMap<usize, Beam>,
+}
+
+impl World {
+    pub fn stiffness_matrix(&self) -> Result<Array<f32, Ix2>, String> {
+        let dof = self.points.len() * 2;
+        let mut a: Vec<f32> = vec![0.0; dof.pow(2)];
+
+        // set up forces exerted by displacements of other beams
+        for (_point_id, point) in &self.points {
+            for beam_id in &point.beams {
+                let beam = self.beams.get(beam_id).unwrap();
+                let stiffness = beam.stiffness(self);
+                let other_point = self.points.get(&beam.other_point(point.id)).unwrap();
+                let ax = 2*(point.id-1);
+                let ay = ax + 1;
+                let bx = 2*(other_point.id-1);
+                let by = bx + 1;
+               
+                let angle = beam.angle(self);
+                println!("point: {} to point {} beam: {} beam_angle: {}", _point_id, other_point.id, beam.id, beam.angle(self).to_degrees());
+                let c2 = angle.cos().powi(2);
+                let s2 = angle.sin().powi(2);
+                let cs = angle.cos() * angle.sin();
+
+                // F_AX += K_AB * c2 * (U_AX)
+                a[ax*dof + ax] += stiffness * c2;
+                // F_AX += K_AB * cs * U_AY
+                a[ax*dof + ay] += stiffness * cs;
+                // F_AY += K_AB * cs * U_AX
+                a[ay*dof + ax] += stiffness * cs;
+                // F_AY += K_AB * s^2 * U_AY
+                a[ay*dof + ay] += stiffness * s2;
+
+                // F_AX += K_AB -c^2 U_BX
+                a[ax*dof + bx] += stiffness * -c2;
+                // F_AX += K_AB -cs U_BY
+                a[ax*dof + by] += stiffness * -cs;
+                // F_AY += K_AB -cs U_BX
+                a[ay*dof + bx] += stiffness * -cs;
+                // F_AY += K_AB -s^2 U_BY
+                a[ay*dof + by] += stiffness * -s2;
+            }
+        }
+
+        match Array2::from_shape_vec((dof, dof), a) {
+            Ok(arr) => Ok(arr),
+            Err(e) => Err(e.to_string()),
+        }
+    }
+
+    pub fn link(&mut self, id1: usize, id2: usize, new_beam: NewBeam) -> Result<(), &str>{
+        let b_id = self.beams.len();
+        let p1 = match self.points.get_mut(&id1) {
+            Some(point) => point,
+            None => return Err("Point with id1 not found"),
+        };
+        p1.beams.push(b_id);
+        let p1_id = p1.id;
+
+        let p2 = match self.points.get_mut(&id2) {
+            Some(point) => point,
+            None => return Err("Point with id1 not found"),
+        };
+        let p2_id = p2.id;
+        p2.beams.push(b_id);
+
+        let mut b = Beam::new(p1_id, p2_id, new_beam);
+        b.id = b_id;
+        self.beams.insert(b_id, b);
+
+        Ok(())
+    }
+}
+
+impl From<Vec<Point>> for World {
+    /// Create a world from a list of points.
+    /// Expects each point to have a unique ID.
+    /// If this is not the case, then some points will be overriden.
+    fn from(points: Vec<Point>) -> World {
+        let mut points = points;
+        let mut w = World { points: HashMap::new(), beams: HashMap::new() };
+
+        while let Some(point) = points.pop() {
+            w.points.insert(point.id, point);
+        }
+
+        w
+    }
+}