Simple truss stiffness matrix calculator in Rust

I recently developed this calculator for some FEA homework i had. I'm still trying to familiarize myself with fundamental rust concepts like the ownership and borrowing. Pls let me know if there is something that could use some work/changes so that I follow rust practices better.

FEA = Finite Element Analysis, the purpose of the code is to produce the stiffness matrix need in the analysis of truss members. Since you have to take into account the angle when calculating the stiffness matrix, theres a lot of tedious computation done by hand involved so this just makes much computation faster.

matrix.rs

use rand::{thread_rng, Rng};

pub struct Matrix {
pub rows: usize,
pub cols: usize,
pub data: Vec<Vec<f64>>
}

impl Matrix {
pub fn zeros(rows: usize, cols: usize) -> Matrix {
Matrix {
rows,
cols,
data: vec![vec![0.0; cols]; rows]
}
}

pub fn value(rows: usize, cols: usize, value: f64) -> Matrix {
Matrix {
rows,
cols,
data: vec![vec![value; cols]; rows]
}
}

pub fn random(rows: usize, cols: usize) -> Matrix {

let mut res = Matrix::zeros(rows, cols);

for i in 0..rows {
for j in 0..cols {
res.data[i][j] = rng.gen::<f64>() * 2.0 - 1.0;
}
}

res
}

pub fn from(data: Vec<Vec<f64>> ) -> Matrix {
Matrix {
rows: data.len(),
cols: data[0].len(),
data,
}
}

pub fn stiffness(modulus: f64, area: f64, length: f64, angle: f64) -> Matrix {
let spring_constant = modulus * area / length;
Matrix::from(vec![
vec![spring_constant * angle.cos().powi(2), spring_constant * angle.cos() * angle.sin(), -spring_constant * angle.cos().powi(2), -spring_constant * angle.cos() * angle.sin()],
vec![spring_constant * angle.cos() * angle.sin(), spring_constant * angle.sin().powi(2), -spring_constant * angle.cos() * angle.sin(), -spring_constant * angle.sin().powi(2)],
vec![-spring_constant * angle.cos().powi(2), -spring_constant * angle.cos() * angle.sin(), spring_constant * angle.cos().powi(2), spring_constant * angle.cos() * angle.sin()],
vec![-spring_constant * angle.cos() * angle.sin(), -spring_constant * angle.sin().powi(2), spring_constant * angle.cos() * angle.sin(), spring_constant * angle.sin().powi(2)],
])
}

pub fn multiply(&mut self, other: &Matrix) -> Matrix {
if self.cols != other.rows {
panic!("Attempted to multiply by matrix of incorrect dimensions.");
}

let mut res = Matrix::zeros(self.rows, other.cols);

for i in 0..self.rows {
for j in 0..other.cols {
let mut sum = 0.0;
for k in 0..self.cols {
sum += self.data[i][k] * other.data[k][j];
}

res.data[i][j] = sum;
}
}

res
}

pub fn add(&mut self, other: &Matrix) -> Matrix {
if self.rows != other.rows || self.cols != other.cols {
panic!("Attempted to add matrix of incorrect dimensions.")
}

let mut res = Matrix::zeros(self.rows, self.cols);

for i in 0..self.rows {
for j in 0..self.cols {
res.data[i][j] = self.data[i][j] + other.data[i][j];
}
}

res
}

pub fn dot_product(&mut self, other: &Matrix) -> Matrix {
if self.rows != other.rows || self.cols != other.cols {
panic!("Attempted to perform dot product on matrix of incorrect dimensions.")
}

let mut res = Matrix::zeros(self.rows, self.cols);

for i in 0..self.rows {
for j in 0..self.cols {
res.data[i][j] = self.data[i][j] * other.data[i][j];
}
}

res
}

pub fn subtract(&mut self, other: &Matrix) -> Matrix {
if self.rows != other.rows || self.cols != other.cols {
panic!("Attempted to subtract matrix of incorrect dimensions.")
}

let mut res = Matrix::zeros(self.rows, self.cols);

for i in 0..self.rows {
for j in 0..self.cols {
res.data[i][j] = self.data[i][j] - other.data[i][j];
}
}

res
}

pub fn map(&mut self, function: &dyn Fn(f64) -> f64) -> Matrix {
Matrix::from(
(self.data)
.clone()
.into_iter()
.map(|row| row.into_iter().map(|value| function(value)).collect())
.collect()
)
}

pub fn transpose(&mut self) -> Matrix {
let mut res = Matrix::zeros(self.cols, self.rows);

for i in 0..self.rows{
for j in 0..self.cols {
res.data[j][i] = self.data[i][j];
}
}

res
}

pub fn to_string(&self) -> String {
let max_width = self.data.iter().flatten().fold(0, |max, val| {
let width = format!("{:.1}", val).len();
max.max(width)
});

let lines: Vec<String> = self.data.iter().map(|row| {
format!(
"| {} |",
row.iter()
.map(|val| format!("{:>width\$.1} ", val, width = max_width))
.collect::<Vec<_>>()
.join(" | ")
)
}).collect();

lines.join("\n")
}

pub fn print(&self) {
println!("{}\n",self.to_string());
}
}


member.rs

use crate::lib::matrix::Matrix;
pub struct TrussMember {
matrix: Matrix,
nodes: Vec<usize>
}

impl TrussMember {
pub fn new(modulus: f64, area: f64, length: f64, angle: f64, nodes: Vec<usize>) -> TrussMember {
TrussMember {
nodes,
}
}

pub fn expand(&self, total_nodes: usize) -> Matrix {
let mut expanded: Matrix = Matrix::zeros(total_nodes * 2, total_nodes * 2);
let mut x = 0;
let mut y = 0;
for i in 0..total_nodes*2 {
for j in 0..total_nodes*2 {
if self.nodes.contains(&(i / 2 + 1)) && self.nodes.contains(&(j / 2 + 1)) {
expanded.data[i][j] = self.matrix.data[x][y];
if y < self.matrix.data.len() {
y += 1;
}
if y + 1 > self.matrix.data[0].len() {
x += 1;
y = 0;
}
}
}
}
expanded
}

pub fn print(&self) {
println!("{}\n", self.matrix.to_string())
}
}


structure.rs

use crate::lib::{member::TrussMember, matrix::Matrix};

pub struct TrussStructure {
members: Vec<TrussMember>,
nodes: Vec<usize>,
}

impl TrussStructure {
pub fn new(members: Vec<TrussMember>, nodes: Vec<usize>) -> TrussStructure {
TrussStructure { members, nodes }
}

pub fn global_stiffness_matrix(&self) -> Matrix {
let total_nodes = self.nodes.len();
let mut gsm = Matrix::zeros(total_nodes * 2, total_nodes * 2);

for member in &self.members {
let mut member_matrix = member.expand(total_nodes);
}

gsm
}

}


main.rs

use lib::member::TrussMember;
use lib::structure::TrussStructure;

fn main() {
let e_mod = 120.0;
let area =  500.0;
let member1 = TrussMember::new(e_mod, area, 1200.0, 180.0, vec![1, 2]);
let member2 = TrussMember::new(e_mod, area, 2000.0, 233.13, vec![1, 3]);
let member3 = TrussMember::new(e_mod, area, 1600.0, 270.0, vec![1, 4]);
let system = TrussStructure::new(vec![member1,member2,member3], vec![1,2,3,4]);
let expanded = system.global_stiffness_matrix();
expanded.print()

}