src/bin/two_d.rs

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

src/two_d.rs

@@ -1,4 +1,14 @@
 pub mod beam_elements;
-pub mod vector;
 
-pub use vector::*;
+#[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

@@ -1,11 +1,25 @@
 mod beam;
-mod boundary_condition;
 mod cross_section;
-mod point;
 mod world;
 
 pub use beam::*;
-pub use boundary_condition::*;
 pub use cross_section::*;
-pub use point::*;
 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

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

src/two_d/beam_elements/boundary_condition.rs

@@ -1,8 +0,0 @@
-use crate::two_d::Vector;
-
-#[derive(Debug)]
-pub enum BoundaryCondition {
-    Free,
-    Fixed,
-    Force(Vector),
-}

src/two_d/beam_elements/cross_section.rs

@@ -9,10 +9,13 @@ pub enum CrossSection {
 
 impl CrossSection {
     pub fn area(&self) -> f32 {
-        match self {
+        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/point.rs

@@ -1,12 +0,0 @@
-use super::super::Vector;
-use super::{Beam, BoundaryCondition};
-use std::cell::RefCell;
-use std::rc::Weak;
-
-#[derive(Debug)]
-pub struct Point {
-    pub id: usize,
-    pub pos: Vector,
-    pub bc: BoundaryCondition,
-    pub beams: Vec<Weak<RefCell<Beam>>>,
-}

src/two_d/beam_elements/world.rs

@@ -1,13 +1,11 @@
 use super::*;
 use ndarray::prelude::*;
-use std::cell::RefCell;
 use std::collections::HashMap;
-use std::rc::Rc;
 
 #[derive(Debug)]
 pub struct World {
-    pub points: HashMap<usize, Rc<RefCell<Point>>>,
+    pub points: HashMap<usize, Point>,
-    pub beams: HashMap<usize, Rc<RefCell<Beam>>>,
+    pub beams: HashMap<usize, Beam>,
 }
 
 impl World {
@@ -15,41 +13,40 @@ impl World {
         let dof = self.points.len() * 2;
         let mut a: Vec<f32> = vec![0.0; dof.pow(2)];
 
-        for point in self.points.values() {
+        // set up forces exerted by displacements of other beams
-            let point = point.borrow();
+        for (_point_id, point) in &self.points {
-            for beam in &point.beams {
+            for beam_id in &point.beams {
-                let beam = beam.upgrade().unwrap();
+                let beam = self.beams.get(beam_id).unwrap();
-                let beam = beam.borrow();
+                let stiffness = beam.stiffness(self);
-                let other_point = beam.other_point(&point).unwrap();
+                let other_point = self.points.get(&beam.other_point(point.id)).unwrap();
-                let other_point = other_point.borrow();
+                let ax = 2*(point.id-1);
-                let ax = 2 * (point.id - 1);
                 let ay = ax + 1;
-                let bx = 2 * (other_point.id - 1);
+                let bx = 2*(other_point.id-1);
                 let by = bx + 1;
                
-                let angle = beam.angle();
+                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();
-                let stiffness = beam.stiffness();
 
                 // F_AX += K_AB * c2 * (U_AX)
-                a[ax * dof + ax] += stiffness * c2;
+                a[ax*dof + ax] += stiffness * c2;
                 // F_AX += K_AB * cs * U_AY
-                a[ax * dof + ay] += stiffness * cs;
+                a[ax*dof + ay] += stiffness * cs;
                 // F_AY += K_AB * cs * U_AX
-                a[ay * dof + ax] += stiffness * cs;
+                a[ay*dof + ax] += stiffness * cs;
                 // F_AY += K_AB * s^2 * U_AY
-                a[ay * dof + ay] += stiffness * s2;
+                a[ay*dof + ay] += stiffness * s2;
 
                 // F_AX += K_AB -c^2 U_BX
-                a[ax * dof + bx] += stiffness * -c2;
+                a[ax*dof + bx] += stiffness * -c2;
                 // F_AX += K_AB -cs U_BY
-                a[ax * dof + by] += stiffness * -cs;
+                a[ax*dof + by] += stiffness * -cs;
                 // F_AY += K_AB -cs U_BX
-                a[ay * dof + bx] += stiffness * -cs;
+                a[ay*dof + bx] += stiffness * -cs;
                 // F_AY += K_AB -s^2 U_BY
-                a[ay * dof + by] += stiffness * -s2;
+                a[ay*dof + by] += stiffness * -s2;
             }
         }
 
@@ -59,30 +56,24 @@ impl World {
         }
     }
 
-    pub fn link(&mut self, id1: usize, id2: usize, new_beam: NewBeam) -> Result<(), &str> {
+    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(&id1) {
+        let p1 = match self.points.get_mut(&id1) {
             Some(point) => point,
             None => return Err("Point with id1 not found"),
         };
-        let p2 = match self.points.get(&id2) {
+        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 b = Rc::new(RefCell::new(Beam::new(
+        let mut b = Beam::new(p1_id, p2_id, new_beam);
-            Rc::downgrade(p1),
+        b.id = b_id;
-            Rc::downgrade(p2),
-            new_beam,
-        )));
-
-        {
-            b.borrow_mut().id = b_id;
-        }
-
-        p1.borrow_mut().beams.push(Rc::downgrade(&b));
-        p2.borrow_mut().beams.push(Rc::downgrade(&b));
-
         self.beams.insert(b_id, b);
 
         Ok(())
@@ -95,13 +86,10 @@ impl From<Vec<Point>> for World {
     /// 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 {
+        let mut w = World { points: HashMap::new(), beams: HashMap::new() };
-            points: HashMap::new(),
-            beams: HashMap::new(),
-        };
 
         while let Some(point) = points.pop() {
-            w.points.insert(point.id, Rc::new(RefCell::new(point)));
+            w.points.insert(point.id, point);
         }
 
         w

src/two_d/vector.rs

@@ -1,12 +0,0 @@
-#[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())
-    }
-}