qtruss 0.9.3

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
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116

//! Describes a truss element.

use maria_linalg::{
    Matrix,
    Vector,
};

use super::PlottableElement;

#[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,
    /// Location of first endpoint
    oneloc: Vector<2>,
    /// Location of second endpoint
    twoloc: Vector<2>,
    /// 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>,
    /// Length (used to determine volume)
    pub length: f64,
}

impl Element {
    /// Constructs a 'zero' element.
    pub fn zero() -> Self {
        Self {
            one: 0,
            two: 0,
            oneloc: Vector::zero(),
            twoloc: Vector::zero(),
            stiffness: Matrix::zero(),
            direction: Vector::zero(),
            length: 0.0,
        }
    }
    
    /// Constructs a new element.
    pub fn new(
        one: usize,
        oneloc: Vector<2>,
        two: usize,
        twoloc: Vector<2>,
    ) -> 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
        // Note that element cross-section and element stiffness are passed at run-time
        
        // 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(1.0 / length);

        Self {
            one,
            two,
            oneloc,
            twoloc,
            stiffness,
            direction,
            length,
        }
    }

    /// Constructs a plottable element.
    pub fn to_plottable(
        &self,
        top: f64,
        left: f64,
        width: f64,
        height: f64,
        stress: f64,
        maxstress: f64,
        area: f64,
        maxarea: f64,
    ) -> PlottableElement {
        PlottableElement {
            x1: ((self.oneloc[0] - left) / width) as f32,
            x2: ((self.oneloc[1] - left) / width) as f32,
            y1: ((height - self.twoloc[0] + top) / height) as f32,
            y2: ((height - self.twoloc[1] + top) / height) as f32,
            stress: (stress / maxstress) as f32,
            thickness: (area / maxarea).sqrt() as f32,
        }
    }
}