From 5b53b52e2fe42c37df8ef397cb46c6fc445bac5f Mon Sep 17 00:00:00 2001
From: Akbar Rahman <hi@alv.cx>
Date: Wed, 27 Dec 2023 17:17:26 +0000
Subject: [PATCH] reduce problem to solve

---
 src/bin/two_d.rs                 |  23 +++++--
 src/lib.rs                       |   6 ++
 src/two_d/beam_elements/world.rs | 102 +++++++++++++++++++++++++++++--
 3 files changed, 121 insertions(+), 10 deletions(-)

diff --git a/src/bin/two_d.rs b/src/bin/two_d.rs
index 08cbe56..e83d2cd 100644
--- a/src/bin/two_d.rs
+++ b/src/bin/two_d.rs
@@ -1,5 +1,6 @@
 use fea::two_d::beam_elements::*;
 use fea::two_d::*;
+use ndarray::prelude::*;
 
 fn main() {
     let student_id_last_digit = 9_f32;
@@ -20,12 +21,11 @@ fn main() {
     let p2 = Point {
         id: 2,
         pos: Vector(0.0, l * theta.sin()),
-        bc: BoundaryCondition::Force(Vector(0.0, 1000.0)),
+        bc: BoundaryCondition::Fixed,
         beams: vec![],
     };
 
     let l_e3 = 2.0 / theta.sin();
-    println!("l_e3 {l_e3}");
 
     let p3 = Point {
         id: 3,
@@ -80,6 +80,21 @@ fn main() {
     )
     .unwrap();
 
-    println!("{:?}", w);
-    println!("{:?}", w.stiffness_matrix().unwrap());
+    println!("w = {:?}", w);
+    println!("w.shape() = {:?}", w.shape());
+    let stiffness = Array2::from_shape_vec(w.shape(), w.stiffness()).unwrap();
+    println!("w.stiffness = \n{:?}", stiffness);
+    println!("w.displacement() = {:?}", w.displacement());
+    println!("w.force() = {:?}", w.force());
+
+    let reduced_stiffness =
+        Array2::from_shape_vec(w.reduced_shape(), w.reduced_stiffness()).unwrap();
+    let reduced_force = Array2::from_shape_vec((1, w.reduced_dof()), w.reduced_force()).unwrap();
+    let reduced_displacement =
+        Array2::from_shape_vec((1, w.reduced_dof()), w.reduced_displacement()).unwrap();
+    println!("w.reduced_stiffness() = \n{:?}", reduced_stiffness);
+    println!("w.reduced_force() = {:?}", reduced_force);
+    println!("w.reduced_displacement() = {:?}", reduced_displacement);
+
+    //println!("soln: {:?}", reduced_force / reduced_stiffness)
 }
diff --git a/src/lib.rs b/src/lib.rs
index 3f2da4b..cfedc29 100644
--- a/src/lib.rs
+++ b/src/lib.rs
@@ -5,3 +5,9 @@ pub struct Material {
     pub youngs_modulus: f32,
     pub yield_stress: f32,
 }
+
+#[derive(Clone, Debug, PartialEq)]
+pub enum Value {
+    Unknown,
+    Known(f32),
+}
diff --git a/src/two_d/beam_elements/world.rs b/src/two_d/beam_elements/world.rs
index c3447b7..9ea20e4 100644
--- a/src/two_d/beam_elements/world.rs
+++ b/src/two_d/beam_elements/world.rs
@@ -1,5 +1,5 @@
 use super::*;
-use ndarray::prelude::*;
+use crate::Value;
 use std::cell::RefCell;
 use std::collections::HashMap;
 use std::rc::Rc;
@@ -11,8 +11,24 @@ pub struct World {
 }
 
 impl World {
-    pub fn stiffness_matrix(&self) -> Result<Array<f32, Ix2>, String> {
-        let dof = self.points.len() * 2;
+    pub fn dof(&self) -> usize {
+        self.points.len() * 2
+    }
+
+    pub fn reduced_dof(&self) -> usize {
+        self.reduced_displacement().len()
+    }
+
+    pub fn shape(&self) -> (usize, usize) {
+        (self.dof(), self.dof())
+    }
+
+    pub fn reduced_shape(&self) -> (usize, usize) {
+        (self.reduced_dof(), self.reduced_dof())
+    }
+    /// Create the stiffness matrix for the system of equations to solve
+    pub fn stiffness(&self) -> Vec<f32> {
+        let dof = self.dof();
         let mut a: Vec<f32> = vec![0.0; dof.pow(2)];
 
         for point in self.points.values() {
@@ -53,10 +69,84 @@ impl World {
             }
         }
 
-        match Array2::from_shape_vec((dof, dof), a) {
-            Ok(arr) => Ok(arr),
-            Err(e) => Err(e.to_string()),
+        a
+    }
+    fn idx_to_coord(&self, idx: usize) -> (usize, usize) {
+        let shape = self.shape();
+
+        let row = idx / shape.0;
+        let col = idx - (row * shape.0);
+
+        (row, col)
+    }
+
+    pub fn reduced_stiffness(&self) -> Vec<f32> {
+        let d = self.displacement();
+
+        self.stiffness()
+            .into_iter()
+            .enumerate()
+            .filter(|&(i, _)| {
+                let (row, col) = self.idx_to_coord(i);
+                //println!("row = {} d[row] = {:?}", row, d[row]);
+                //println!("col = {} d[col] = {:?}", col, d[col]);
+                d[row] == Value::Unknown && d[col] == Value::Unknown
+            })
+            .map(|(_, e)| e)
+            .collect()
+    }
+
+    pub fn displacement(&self) -> Vec<Value> {
+        let mut a = vec![Value::Known(0_f32); self.dof()];
+
+        for point in self.points.values() {
+            let point = point.borrow();
+            let ax = 2 * (point.id - 1);
+            let ay = ax + 1;
+
+            if let BoundaryCondition::Fixed = point.bc {
+                continue;
+            }
+
+            a[ax] = Value::Unknown;
+            a[ay] = Value::Unknown;
         }
+
+        a
+    }
+
+    pub fn reduced_displacement(&self) -> Vec<Value> {
+        self.displacement()
+            .into_iter()
+            .filter(|d| *d == Value::Unknown)
+            .collect::<Vec<_>>()
+    }
+
+    pub fn force(&self) -> Vec<Value> {
+        let mut a = vec![Value::Unknown; self.dof()];
+        for point in self.points.values() {
+            let point = point.borrow();
+            let ax = 2 * (point.id - 1);
+            let ay = ax + 1;
+
+            if let BoundaryCondition::Force(v) = &point.bc {
+                a[ax] = Value::Known(v.0);
+                a[ay] = Value::Known(v.1);
+            }
+        }
+
+        a
+    }
+
+    pub fn reduced_force(&self) -> Vec<Value> {
+        let d = self.displacement();
+
+        self.force()
+            .into_iter()
+            .enumerate()
+            .filter(|&(i, _)| d[i] == Value::Unknown)
+            .map(|(_, f)| f)
+            .collect()
     }
 
     pub fn link(&mut self, id1: usize, id2: usize, new_beam: NewBeam) -> Result<(), &str> {