This module defines the trait FiniteElement, and provides methods to solve finite element systems.

Example

use finiteelement::*;
//Define a struct that implements the FiniteElement trait


#[derive(Clone)]
pub struct Spring {
   pub a: usize,
   pub b: usize,
   pub l: f64,
   pub k: f64,
}

impl FiniteElement<f64> for Spring {
   fn forces(&self, positions: &[Point<f64>], forces: &mut [f64]) {
       // add to both a and b the force resulting from this spring.
       let ab = positions[self.b].clone() - positions[self.a].clone();
       let norm = ab.norm();
       forces[3 * self.a] += self.k * (norm - self.l) * ab.x / norm;
       forces[3 * self.a + 1] += self.k * (norm - self.l) * ab.y / norm;
       forces[3 * self.a + 2] += self.k * (norm - self.l) * ab.z / norm;
       forces[3 * self.b] -= self.k * (norm - self.l) * ab.x / norm;
       forces[3 * self.b + 1] -= self.k * (norm - self.l) * ab.y / norm;
       forces[3 * self.b + 2] -= self.k * (norm - self.l) * ab.z / norm;
   }

   fn jacobian(&self, positions: &[Point<f64>], jacobian: &mut Sparse<f64>) {
       // add to both a and b the force resulting from this self.
       let ab = positions[self.b].clone() - positions[self.a].clone();
       let norm = ab.norm();
       let norm3 = norm * norm * norm;
       for u in 0..3 {
           for v in 0..3 {
               let j = if u == v {
                   self.k * (1. - self.l / norm + self.l * ab[u] * ab[u] / norm3)
               } else {
                   self.k * self.l * ab[u] * ab[v] / norm3
               };
               // Change in the force on B when moving B.
               jacobian[3 * self.b + u][3 * self.b + v] -= j;
               // Change in the force on A when moving B.
               jacobian[3 * self.a + u][3 * self.b + v] += j;
               // Change in the force on B when moving A.
               jacobian[3 * self.b + u][3 * self.a + v] += j;
               // Change in the force on A when moving A.
               jacobian[3 * self.a + u][3 * self.a + v] -= j;
           }
       }
   }
}

let elts = [
   Spring {
       a: 0,
       b: 1,
       l: 1.,
       k: 1.,
   },
   Spring {
       a: 1,
       b: 2,
       l: 2.,
       k: 0.5,
   },
   Spring {
       a: 0,
       b: 2,
       l: 3.,
       k: 5.
   },
];
let system = (0..elts.len()).map(|i| {Box::new(elts[i].clone()) as Box<dyn FiniteElement<f64>>}).collect::<Vec<_>>();
let mut positions = vec![
 Point { x: 0., y: 0., z: 0. },
 Point {x : 1., y: 0., z: 0.},
 Point{x: 0., y: 1., z: 1.}];

let mut ws = FesWorkspace::new(positions.len());
let epsilon_stop = 1e-4;
let gradient_switch = 1e-3;
let mut solved = false;
for i in (0..20) {
  solved = fes_one_step(&system, &mut positions, epsilon_stop, gradient_switch, &mut ws);
  if solved {
     break;
  }
}
assert!(solved);

Use of the macro provided by finiteelement_macro

To use one of the macro provided by finiteelement_marco, call it without parameter at the begining or your code. The piece of code surrounded by //======== in the example above can be replaced by a call to finiteelement_macros::auto_impl_spring!{}

 use finiteelement::*;
 use std::borrow::Borrow;
auto_impl_spring!{}
pub fn main() {
 let elts = [
     Spring {
         a: 0,
         b: 1,
         l: 1.,
         k: 1.,
     },
     Spring {
         a: 1,
         b: 2,
         l: 2.,
         k: 0.5,
     },
     Spring {
         a: 0,
         b: 2,
         l: 3.,
         k: 5.
     },
 ];
 let system = (0..elts.len()).map(|i| {Box::new(elts[i].clone()) as Box<dyn FiniteElement<f64>>}).collect::<Vec<_>>();
 let mut positions = vec![
   Point { x: 0., y: 0., z: 0. },
   Point {x : 1., y: 0., z: 0.},
   Point{x: 0., y: 1., z: 1.}];

 let mut ws = FesWorkspace::new(positions.len());
 let epsilon_stop = 1e-4;
 let gradient_switch = 1e-3;
 let mut solved = false;
 for i in (0..20) {
    solved = fes_one_step(&system, &mut positions, epsilon_stop, gradient_switch, &mut ws);
    if solved {
       break;
    }
 }
assert!(solved);
}