qtruss 0.2.0

A simple finite-element solver for trusses.
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89

//! Describes a truss element.

use super::{
    Matrix,
    Vector,
};

#[derive(Clone, Copy, Debug)]
/// A single 2D truss element.
pub struct Element {
    /// Node number of first endpoint
    pub one: usize,
    /// Node number of second endpoint
    pub two: usize,
    /// Local stiffness matrix, expressed in global coordinates
    pub stiffness: Matrix<4>,
    /// Unit vector indicating element direction (used to determine if element is in tension or compression)
    pub direction: Vector<2>,
    /// Cross sectional area (used to determine stress)
    pub area: f64,
    /// Length (used to determine volume)
    pub length: f64,
}

impl Element {
    /// Constructs a 'zero' element.
    pub fn zero() -> Self {
        Self {
            one: 0,
            two: 0,
            stiffness: Matrix::zero(),
            direction: Vector::zero(),
            area: 0.0,
            length: 0.0,
        }
    }
    
    /// Constructs a new element.
    pub fn new(
        one: usize,
        oneloc: Vector<2>,
        two: usize,
        twoloc: Vector<2>,
        e: f64,
        area: f64,
    ) -> Self {
        // Compute element length
        let r = twoloc - oneloc;
        let length = r.norm();
        let direction = r.scale(1.0 / length);
        
        // Compute element angle
        let ex = Vector::new([1.0, 0.0]);
        let ey = Vector::new([0.0, 1.0]);
        let opposite = r.dot(ey);
        let adjacent = r.dot(ex);
        let hypotenuse = (opposite.powf(2.0) + adjacent.powf(2.0)).sqrt();
        let c = adjacent / hypotenuse;
        let s = opposite / hypotenuse;

        // Compute element stiffness matrix in global coordinates
        
        // The Stiffness Matrix
        // (Multiply against displacement vector to get force vector)
        // 
        //    1X 1Y 2X 2Y
        // 1X  x  x  x  x
        // 1Y  x  x  x  x
        // 2X  x  x  x  x
        // 2Y  x  x  x  x
        // 
        // Element `(i, j)` is the coefficient of displacement `j` on force `i`
        let stiffness = Matrix::new([
            [  c * c,    c * s,  - c * c,  - c * s],
            [  c * s,    s * s,  - c * s,  - s * s],
            [- c * c,  - c * s,    c * c,    c * s],
            [- c * s,  - s * s,    c * s,    s * s],
        ]).scale(e * area / length);

        Self {
            one,
            two,
            stiffness,
            direction,
            area,
            length,
        }
    }
}